

<!DOCTYPE html>
<!--[if IE 8]><html class="no-js lt-ie9" lang="en" > <![endif]-->
<!--[if gt IE 8]><!--> <html class="no-js" lang="en" > <!--<![endif]-->
<head>
  <meta charset="utf-8">
  
  <meta name="viewport" content="width=device-width, initial-scale=1.0">
  
  <title>quippy.crack &mdash; quippy 7ea04c5+ documentation</title>
  

  
  
    <link rel="shortcut icon" href="../../_static/favicon.ico"/>
  
  
  

  
  <script type="text/javascript" src="../../_static/js/modernizr.min.js"></script>
  
    
      <script type="text/javascript" id="documentation_options" data-url_root="../../" src="../../_static/documentation_options.js"></script>
        <script type="text/javascript" src="../../_static/jquery.js"></script>
        <script type="text/javascript" src="../../_static/underscore.js"></script>
        <script type="text/javascript" src="../../_static/doctools.js"></script>
        <script type="text/javascript" src="../../_static/language_data.js"></script>
        <script async="async" type="text/javascript" src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/latest.js?config=TeX-AMS-MML_HTMLorMML"></script>
        <script type="text/x-mathjax-config">MathJax.Hub.Config({"tex2jax": {"processClass": "math|output_area", "processEscapes": true, "ignoreClass": "document", "inlineMath": [["$", "$"], ["\\(", "\\)"]]}})</script>
    
    <script type="text/javascript" src="../../_static/js/theme.js"></script>

    

  
  <link rel="stylesheet" href="../../_static/css/theme.css" type="text/css" />
  <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
  <link rel="stylesheet" href="../../_static/theme_overrides.css" type="text/css" />
    <link rel="index" title="Index" href="../../genindex.html" />
    <link rel="search" title="Search" href="../../search.html" /> 
</head>

<body class="wy-body-for-nav">

   
  <div class="wy-grid-for-nav">
    
    <nav data-toggle="wy-nav-shift" class="wy-nav-side">
      <div class="wy-side-scroll">
        <div class="wy-side-nav-search" >
          

          
            <a href="../../index.html" class="icon icon-home"> quippy
          

          
            
            <img src="../../_static/hybrid.png" class="logo" alt="Logo"/>
          
          </a>

          
            
            
              <div class="version">
                7ea04c5+
              </div>
            
          

          
<div role="search">
  <form id="rtd-search-form" class="wy-form" action="../../search.html" method="get">
    <input type="text" name="q" placeholder="Search docs" />
    <input type="hidden" name="check_keywords" value="yes" />
    <input type="hidden" name="area" value="default" />
  </form>
</div>

          
        </div>

        <div class="wy-menu wy-menu-vertical" data-spy="affix" role="navigation" aria-label="main navigation">
          
    
            
            
              
            
            
              <ul>
<li class="toctree-l1"><a class="reference internal" href="../../intro.html">Introduction to <code class="docutils literal notranslate"><span class="pre">quippy</span></code></a></li>
<li class="toctree-l1"><a class="reference internal" href="../../install.html">Installation of QUIP and quippy</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../Tutorials/index.html">Tutorials</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../quippy.html">Quippy library reference</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../visualisation.html">Visualisation Tools</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../fortran_wrapper.html">Appendix: wrapping Fortran 90 Code</a></li>
<li class="toctree-l1"><a class="reference internal" href="../../teach_sparse.html">Command line options of the teach_sparse main program</a></li>
</ul>

            
          
    <p></p>
    <ul>
        <li><a href="../../genindex.html">Index</a></li>
        <li><a href="../../py-modindex.html">Module Index</a></li>
        <li><a href="../../search.html">Search Page</a></li>
        <li><a href="../index.html">Module Source Code</a></li>
    </ul>
  
        </div>
      </div>
    </nav>

    <section data-toggle="wy-nav-shift" class="wy-nav-content-wrap">

      
      <nav class="wy-nav-top" aria-label="top navigation">
        
          <i data-toggle="wy-nav-top" class="fa fa-bars"></i>
          <a href="../../index.html">quippy</a>
        
      </nav>


      <div class="wy-nav-content">
        
        <div class="rst-content">
        
          















<div role="navigation" aria-label="breadcrumbs navigation">

  <ul class="wy-breadcrumbs">
    
      <li><a href="../../index.html">Docs</a> &raquo;</li>
        
          <li><a href="../index.html">Module code</a> &raquo;</li>
        
      <li>quippy.crack</li>
    
    
      <li class="wy-breadcrumbs-aside">
        
      </li>
    
  </ul>

  
  <hr/>
</div>
          <div role="main" class="document" itemscope="itemscope" itemtype="http://schema.org/Article">
           <div itemprop="articleBody">
            
  <h1>Source code for quippy.crack</h1><div class="highlight"><pre>
<span></span><span class="c1"># HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   quippy: Python interface to QUIP atomistic simulation library</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   Copyright James Kermode 2010</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   These portions of the source code are released under the GNU General</span>
<span class="c1"># HQ X   Public License, version 2, http://www.gnu.org/copyleft/gpl.html</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   If you would like to license the source code under different terms,</span>
<span class="c1"># HQ X   please contact James Kermode, james.kermode@gmail.com</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   When using this software, please cite the following reference:</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ X   http://www.jrkermode.co.uk/quippy</span>
<span class="c1"># HQ X</span>
<span class="c1"># HQ XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX</span>

<span class="kn">import</span> <span class="nn">os</span><span class="o">,</span> <span class="nn">warnings</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="k">import</span> <span class="n">sqrt</span><span class="p">,</span> <span class="n">pi</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">numpy.ma</span> <span class="k">as</span> <span class="nn">ma</span>
<span class="kn">from</span> <span class="nn">quippy.mpi_context</span> <span class="k">import</span> <span class="n">MPI_context</span>
<span class="kn">from</span> <span class="nn">quippy.potential</span> <span class="k">import</span> <span class="n">Potential</span>
<span class="kn">from</span> <span class="nn">quippy.atoms</span> <span class="k">import</span> <span class="n">Atoms</span>
<span class="kn">from</span> <span class="nn">quippy.structures</span> <span class="k">import</span> <span class="n">transform</span><span class="p">,</span> <span class="n">supercell</span><span class="p">,</span> <span class="n">MillerPlane</span><span class="p">,</span> <span class="n">MillerDirection</span>
<span class="kn">from</span> <span class="nn">quippy.system</span> <span class="k">import</span> <span class="n">print_title</span><span class="p">,</span> <span class="n">verbosity_push</span>
<span class="kn">from</span> <span class="nn">quippy.clusters</span> <span class="k">import</span> <span class="n">HYBRID_NO_MARK</span>
<span class="kn">from</span> <span class="nn">quippy.units</span> <span class="k">import</span> <span class="n">EV_A3_IN_GPA</span>
<span class="kn">from</span> <span class="nn">quippy.surface</span> <span class="k">import</span> <span class="n">J_PER_M2</span>
<span class="kn">from</span> <span class="nn">quippy.farray</span> <span class="k">import</span> <span class="n">fzeros</span><span class="p">,</span> <span class="n">farray</span>

<span class="n">MPA_SQRT_M</span> <span class="o">=</span> <span class="mf">1e-3</span><span class="o">/</span><span class="n">EV_A3_IN_GPA</span><span class="o">*</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">1.0e10</span><span class="p">)</span>

<span class="n">__all__</span> <span class="o">=</span> <span class="p">[]</span>

<span class="k">try</span><span class="p">:</span>
   <span class="kn">import</span> <span class="nn">quippy.cracktools</span><span class="o">,</span> <span class="nn">quippy.crackparams</span>
   <span class="kn">from</span> <span class="nn">quippy.crackparams</span> <span class="k">import</span> <span class="o">*</span>
   <span class="kn">from</span> <span class="nn">quippy.cracktools</span> <span class="k">import</span> <span class="o">*</span>
   <span class="n">__all__</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">quippy</span><span class="o">.</span><span class="n">crackparams</span><span class="o">.</span><span class="n">__all__</span><span class="p">)</span>
   <span class="n">__all__</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="n">quippy</span><span class="o">.</span><span class="n">cracktools</span><span class="o">.</span><span class="n">__all__</span><span class="p">)</span>

<span class="k">except</span> <span class="ne">ImportError</span><span class="p">:</span>
   <span class="n">warnings</span><span class="o">.</span><span class="n">warn</span><span class="p">(</span><span class="s1">&#39;crack utilities not available&#39;</span><span class="p">)</span>

<span class="n">__all__</span><span class="o">.</span><span class="n">extend</span><span class="p">([</span><span class="s1">&#39;crack_rescale_homogeneous_xy&#39;</span><span class="p">,</span>
                <span class="s1">&#39;crack_rescale_uniaxial&#39;</span><span class="p">,</span>
                <span class="s1">&#39;makecrack&#39;</span><span class="p">,</span>
                <span class="s1">&#39;crack_strain_energy_release_rate&#39;</span><span class="p">,</span>
                <span class="s1">&#39;crack_strain&#39;</span><span class="p">,</span>
                <span class="s1">&#39;crack_find_griffith_load&#39;</span><span class="p">,</span>
                <span class="s1">&#39;stress_intensity_factor&#39;</span><span class="p">,</span>
                <span class="s1">&#39;make_crack_advance_map&#39;</span><span class="p">,</span>
                <span class="s1">&#39;find_crack_tip_coordination&#39;</span><span class="p">,</span>
                <span class="s1">&#39;irwin_modeI_crack_tip_stress_field&#39;</span><span class="p">,</span>
                <span class="s1">&#39;strain_to_G&#39;</span><span class="p">,</span>
                <span class="s1">&#39;G_to_strain&#39;</span><span class="p">,</span>
                <span class="s1">&#39;get_strain&#39;</span><span class="p">,</span>
                <span class="s1">&#39;get_energy_release_rate&#39;</span><span class="p">,</span>
                <span class="s1">&#39;get_stress_intensity_factor&#39;</span><span class="p">,</span>
                <span class="s1">&#39;fit_crack_stress_field&#39;</span><span class="p">,</span>
                <span class="s1">&#39;find_crack_tip_stress_field&#39;</span><span class="p">,</span>
                <span class="s1">&#39;plot_stress_fields&#39;</span><span class="p">,</span>
                <span class="s1">&#39;thin_strip_displacement_y&#39;</span><span class="p">,</span>
                <span class="s1">&#39;print_crack_system&#39;</span><span class="p">,</span>
                <span class="s1">&#39;ConstantStrainRate&#39;</span><span class="p">])</span>

<span class="k">def</span> <span class="nf">makecrack_main</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;Given a CrackParams object `param`, construct and return a new crack slab Atoms object.&quot;&quot;&quot;</span>

    <span class="n">xmlfilename</span> <span class="o">=</span> <span class="n">stem</span><span class="o">+</span><span class="s1">&#39;.xml&#39;</span>

    <span class="n">print_title</span><span class="p">(</span><span class="s1">&#39;Initialisation&#39;</span><span class="p">)</span>

    <span class="n">verbosity_push</span><span class="p">(</span><span class="n">params</span><span class="o">.</span><span class="n">io_verbosity</span><span class="p">)</span>
    <span class="n">params</span><span class="o">.</span><span class="n">print_</span><span class="p">()</span>

    <span class="nb">print</span><span class="p">(</span><span class="s2">&quot;Initialising classical potential with args &quot;</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">classical_args</span><span class="o">.</span><span class="n">strip</span><span class="p">()</span> <span class="o">+</span>
          <span class="s2">&quot; from file &quot;</span> <span class="o">+</span> <span class="n">xmlfilename</span><span class="p">)</span>
    <span class="n">classicalpot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">params</span><span class="o">.</span><span class="n">classical_args</span><span class="p">,</span> <span class="n">param_filename</span><span class="o">=</span><span class="n">xmlfilename</span><span class="p">)</span>
    <span class="n">classicalpot</span><span class="o">.</span><span class="n">print_</span><span class="p">()</span>

    <span class="n">mpi_glob</span> <span class="o">=</span> <span class="n">MPI_context</span><span class="p">()</span>

    <span class="n">crack_slab</span><span class="p">,</span> <span class="n">width</span><span class="p">,</span> <span class="n">height</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="n">v</span><span class="p">,</span> <span class="n">v2</span><span class="p">,</span> <span class="n">bulk</span> <span class="o">=</span> <span class="n">crack_make_slab</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">classicalpot</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_free_surfaces</span><span class="p">:</span>
       <span class="n">depth</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="k">else</span><span class="p">:</span>
       <span class="n">depth</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span>

    <span class="c1"># Save bulk cube (used for qm_rescale_r parameter in crack code)</span>
    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">qm_args</span><span class="o">.</span><span class="n">startswith</span><span class="p">(</span><span class="s1">&#39;TB&#39;</span><span class="p">):</span>
        <span class="n">bigger_bulk</span> <span class="o">=</span> <span class="n">supercell</span><span class="p">(</span><span class="n">bulk</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">2</span><span class="p">)</span>
        <span class="n">bulk</span> <span class="o">=</span> <span class="n">bigger_bulk</span>
    <span class="n">bulk</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">stem</span><span class="o">+</span><span class="s1">&#39;_bulk.xyz&#39;</span><span class="p">)</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">stem</span><span class="o">+</span><span class="s1">&#39;_slab.xyz&#39;</span><span class="p">)</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigWidth&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">width</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">height</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigDepth&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">depth</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;YoungsModulus&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">E</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yx&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">v</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yz&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">v2</span>

    <span class="c1"># Open surfaces, remain periodic in z direction (normal to plane)</span>
    <span class="c1"># and optionally also in x direction if crack_double_ended is true</span>
    <span class="k">if</span> <span class="ow">not</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_double_ended</span><span class="p">:</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">set_lattice</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">,</span> <span class="kc">False</span><span class="p">)</span>

    <span class="c1"># 3D crack with free surfaces at z = +/- depth/2</span>
    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_free_surfaces</span><span class="p">:</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span> <span class="o">-=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">.</span><span class="n">mean</span><span class="p">()</span> <span class="c1"># center on z=0</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">set_lattice</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">,</span> <span class="kc">False</span><span class="p">)</span>

    <span class="c1"># Map atoms into cell AFTER changing to the new lattice</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">map_into_cell</span><span class="p">()</span>

    <span class="n">miny</span><span class="p">,</span> <span class="n">maxy</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
    <span class="k">assert</span> <span class="nb">abs</span><span class="p">((</span><span class="n">maxy</span><span class="o">-</span><span class="n">miny</span><span class="p">)</span> <span class="o">-</span> <span class="n">height</span><span class="p">)</span> <span class="o">&lt;</span> <span class="mf">1e-5</span>  <span class="c1"># be sure that remapping didn&#39;t change height of slab</span>

    <span class="c1"># Add various properties to crack_slab</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;hybrid&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;hybrid_mark&#39;</span><span class="p">,</span> <span class="n">HYBRID_NO_MARK</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;changed_nn&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;move_mask&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;nn&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;old_nn&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;md_old_changed_nn&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;edge_mask&#39;</span><span class="p">,</span> <span class="mi">0</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;crack_surface&#39;</span><span class="p">,</span> <span class="kc">False</span><span class="p">)</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;crack_front&#39;</span><span class="p">,</span> <span class="kc">False</span><span class="p">)</span>
    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles</span><span class="p">:</span>
	<span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;fixdip&#39;</span><span class="p">,</span> <span class="kc">False</span><span class="p">)</span>

    <span class="n">print_title</span><span class="p">(</span><span class="s1">&#39;Fixing Atoms&#39;</span><span class="p">)</span>

    <span class="c1"># Fix top and bottom edges - anything within crack_edge_fix_tol of ymax or ymin is fixed</span>

    <span class="n">miny</span><span class="p">,</span> <span class="n">maxy</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">move_mask</span><span class="p">[:]</span> <span class="o">=</span> <span class="mi">1</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">move_mask</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxy</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_edge_fix_tol</span><span class="p">)</span> <span class="o">|</span>
                         <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">miny</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_edge_fix_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>

    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_sides</span><span class="p">:</span>
        <span class="n">maxx</span><span class="p">,</span> <span class="n">minx</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">move_mask</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_edge_fix_tol</span><span class="p">)</span> <span class="o">|</span>
                             <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">minx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_edge_fix_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">0</span>


    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;</span><span class="si">%d</span><span class="s1"> atoms. </span><span class="si">%d</span><span class="s1"> fixed atoms&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">n</span><span class="p">,</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">n</span> <span class="o">-</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">move_mask</span><span class="o">.</span><span class="n">sum</span><span class="p">()))</span>

    <span class="n">print_title</span><span class="p">(</span><span class="s1">&#39;Setting edge mask&#39;</span><span class="p">)</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">edge_mask</span><span class="p">[:]</span> <span class="o">=</span> <span class="mi">0</span>

    <span class="n">minx</span><span class="p">,</span> <span class="n">maxx</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">edge_mask</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">minx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)</span> <span class="o">|</span>
                         <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">1</span>

    <span class="n">miny</span><span class="p">,</span> <span class="n">maxy</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
    <span class="n">crack_slab</span><span class="o">.</span><span class="n">edge_mask</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">miny</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)</span> <span class="o">|</span>
                         <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxy</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">1</span>

    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_free_surfaces</span><span class="p">:</span>
        <span class="c1"># Open surfaces at +/- z</span>
        <span class="n">minz</span><span class="p">,</span> <span class="n">maxz</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">edge_mask</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">-</span><span class="n">minz</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)</span> <span class="o">|</span>
                             <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxz</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">selection_edge_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">1</span>

    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles</span><span class="p">:</span>
        <span class="n">print_title</span><span class="p">(</span><span class="s1">&#39;Fixing dipoles&#39;</span><span class="p">)</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">fixdip</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxy</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles_tol</span><span class="p">)</span> <span class="o">|</span>
                          <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">miny</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">1</span>
        <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_sides</span><span class="p">:</span>
                <span class="n">maxx</span><span class="p">,</span> <span class="n">minx</span> <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
                <span class="n">crack_slab</span><span class="o">.</span><span class="n">fixdip</span><span class="p">[(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">maxx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles_tol</span><span class="p">)</span> <span class="o">|</span>
                                  <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span><span class="o">-</span><span class="n">minx</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_fix_dipoles_tol</span><span class="p">)]</span> <span class="o">=</span> <span class="mi">1</span>


    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_curved_front</span><span class="p">:</span>
       <span class="n">crack_make_seed_curved_front</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
       <span class="n">crack_make_seed</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>
       <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_apply_initial_load</span><span class="p">:</span>
          <span class="n">crack_calc_load_field</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">,</span> <span class="n">classicalpot</span><span class="p">,</span>
                              <span class="n">params</span><span class="o">.</span><span class="n">crack_loading</span><span class="p">,</span> <span class="n">overwrite_pos</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
                              <span class="n">mpi</span><span class="o">=</span><span class="n">mpi_glob</span><span class="p">)</span>

    <span class="n">crack_slab</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="s1">&#39;dump.xyz&#39;</span><span class="p">)</span>
    <span class="n">crack_update_connect</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>

    <span class="k">if</span> <span class="ow">not</span> <span class="n">params</span><span class="o">.</span><span class="n">simulation_classical</span><span class="p">:</span>
        <span class="k">if</span> <span class="p">(</span><span class="n">params</span><span class="o">.</span><span class="n">selection_method</span><span class="o">.</span><span class="n">strip</span><span class="p">()</span> <span class="o">==</span> <span class="s1">&#39;crack_front&#39;</span> <span class="ow">or</span>
            <span class="n">params</span><span class="o">.</span><span class="n">crack_tip_method</span><span class="o">.</span><span class="n">strip</span><span class="p">()</span> <span class="o">==</span> <span class="s1">&#39;local_energy&#39;</span><span class="p">):</span>
            <span class="n">classicalpot</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">local_energy</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>

        <span class="n">crack_setup_marks</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>
        <span class="n">crack_update_selection</span><span class="p">(</span><span class="n">crack_slab</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">any_per_atom_tau</span><span class="p">():</span>
        <span class="c1"># Set up per_atom_tau property for ramped Langevin thermostat:</span>
        <span class="c1">#</span>
        <span class="c1">#    tau</span>
        <span class="c1">#    ^</span>
        <span class="c1">#    |\        /|                     |\        /|  max_tau</span>
        <span class="c1">#    | \      / |                     | \      / |</span>
        <span class="c1">#    |  \    /  |     constant E      |  \    /  |</span>
        <span class="c1">#    |   \  /   |      (tau = 0)      |   \  /   |</span>
        <span class="c1">#    |    \/    |                     |    \/    |</span>
        <span class="c1">#    +----------+---------------------+----------+---&gt; x</span>
        <span class="c1">#   -w/2      -w/2+r                 w/2-r      w/2</span>

        <span class="n">w_by_2</span>   <span class="o">=</span> <span class="n">crack_slab</span><span class="o">.</span><span class="n">OrigWidth</span><span class="o">/</span><span class="mf">2.</span>
        <span class="n">ramp_len</span> <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_thermostat_ramp_length</span>
        <span class="n">max_tau</span>  <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_thermostat_ramp_max_tau</span>
        <span class="nb">print</span> <span class="s1">&#39;Adding thermostat ramp with length&#39;</span><span class="p">,</span> <span class="n">ramp_len</span><span class="p">,</span> <span class="s1">&#39;max_tau&#39;</span><span class="p">,</span> <span class="n">max_tau</span>

        <span class="nd">@np</span><span class="o">.</span><span class="n">vectorize</span>
        <span class="k">def</span> <span class="nf">tau</span><span class="p">(</span><span class="n">x</span><span class="p">):</span>
            <span class="k">if</span> <span class="n">x</span> <span class="o">&lt;</span> <span class="o">-</span><span class="n">w_by_2</span> <span class="o">+</span> <span class="n">ramp_len</span><span class="o">/</span><span class="mi">2</span><span class="p">:</span>
                <span class="n">q</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span><span class="o">+</span><span class="n">w_by_2</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
                <span class="k">return</span> <span class="n">max_tau</span><span class="o">*</span><span class="p">(</span><span class="mf">1.</span><span class="o">-</span> <span class="n">q</span><span class="p">)</span>
            <span class="k">elif</span> <span class="p">(</span><span class="n">x</span> <span class="o">&gt;</span> <span class="o">-</span><span class="n">w_by_2</span> <span class="o">+</span> <span class="n">ramp_len</span><span class="o">/</span><span class="mi">2</span> <span class="ow">and</span>
                  <span class="n">x</span> <span class="o">&lt;</span> <span class="o">-</span><span class="n">w_by_2</span> <span class="o">+</span> <span class="n">ramp_len</span><span class="p">):</span>
                <span class="n">q</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span><span class="o">+</span><span class="n">w_by_2</span><span class="o">-</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
                <span class="k">return</span> <span class="n">max_tau</span><span class="o">*</span><span class="n">q</span>
            <span class="k">elif</span> <span class="p">(</span><span class="n">x</span> <span class="o">&gt;</span> <span class="o">-</span><span class="n">w_by_2</span> <span class="o">+</span> <span class="n">ramp_len</span> <span class="ow">and</span>
                  <span class="n">x</span> <span class="o">&lt;</span> <span class="n">w_by_2</span> <span class="o">-</span> <span class="n">ramp_len</span><span class="p">):</span>
                <span class="k">return</span> <span class="mf">0.</span>
            <span class="k">elif</span> <span class="p">(</span><span class="n">x</span> <span class="o">&gt;</span> <span class="n">w_by_2</span> <span class="o">-</span> <span class="n">ramp_len</span> <span class="ow">and</span>
                  <span class="n">x</span> <span class="o">&lt;</span> <span class="n">w_by_2</span> <span class="o">-</span> <span class="n">ramp_len</span><span class="o">/</span><span class="mi">2</span><span class="p">):</span>
                <span class="n">q</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">w_by_2</span><span class="o">+</span><span class="n">ramp_len</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
                <span class="k">return</span> <span class="n">max_tau</span><span class="o">*</span><span class="p">(</span><span class="mf">1.</span><span class="o">-</span> <span class="n">q</span><span class="p">)</span>
            <span class="k">else</span><span class="p">:</span>
                <span class="n">q</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">w_by_2</span><span class="o">+</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">ramp_len</span><span class="o">/</span><span class="mf">2.</span><span class="p">)</span>
                <span class="k">return</span> <span class="n">max_tau</span><span class="o">*</span><span class="n">q</span>
        <span class="n">crack_slab</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;per_atom_tau&#39;</span><span class="p">,</span> <span class="n">tau</span><span class="p">(</span><span class="n">crack_slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]))</span>

    <span class="k">return</span> <span class="n">crack_slab</span>


<span class="k">try</span><span class="p">:</span>
   <span class="kn">from</span> <span class="nn">quippy.cracktools</span> <span class="k">import</span> <span class="n">crack_measure_g</span> <span class="k">as</span> <span class="n">_crack_measure_g</span>

   <span class="k">def</span> <span class="nf">crack_measure_g</span><span class="p">(</span><span class="n">crack</span><span class="p">,</span> <span class="n">YoungsModulus</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">PoissonRatio_yx</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">OrigHeight</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
       <span class="k">if</span> <span class="n">YoungsModulus</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
           <span class="n">YoungsModulus</span> <span class="o">=</span> <span class="n">crack</span><span class="o">.</span><span class="n">YoungsModulus</span>
       <span class="k">if</span> <span class="n">PoissonRatio_yx</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
           <span class="n">PoissonRatio_yx</span> <span class="o">=</span> <span class="n">crack</span><span class="o">.</span><span class="n">PoissonRatio_yx</span>
       <span class="k">if</span> <span class="n">OrigHeight</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
           <span class="n">OrigHeight</span> <span class="o">=</span> <span class="n">crack</span><span class="o">.</span><span class="n">OrigHeight</span>
       <span class="k">return</span> <span class="n">_crack_measure_g</span><span class="p">(</span><span class="n">crack</span><span class="p">,</span> <span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">OrigHeight</span><span class="p">)</span>

   <span class="kn">import</span> <span class="nn">functools</span>
   <span class="n">crack_measure_g</span> <span class="o">=</span> <span class="n">functools</span><span class="o">.</span><span class="n">update_wrapper</span><span class="p">(</span><span class="n">crack_measure_g</span><span class="p">,</span> <span class="n">quippy</span><span class="o">.</span><span class="n">cracktools</span><span class="o">.</span><span class="n">crack_measure_g</span><span class="p">)</span>

<span class="k">except</span> <span class="ne">ImportError</span><span class="p">:</span>
   <span class="k">pass</span>

<div class="viewcode-block" id="stress_intensity_factor"><a class="viewcode-back" href="../../crack.html#quippy.crack.stress_intensity_factor">[docs]</a><span class="k">def</span> <span class="nf">stress_intensity_factor</span><span class="p">(</span><span class="n">at</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Returns stress instensity factor for mode I loading (K_I) in MPa \sqrt(m)</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">G</span> <span class="o">=</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">at</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">crack_g_to_k</span><span class="p">(</span><span class="n">G</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="s1">&#39;plane strain&#39;</span><span class="p">)</span><span class="o">/</span><span class="mf">1.0e6</span></div>

<span class="k">def</span> <span class="nf">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">params</span><span class="p">,</span> <span class="n">strain</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">G</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
    <span class="k">if</span> <span class="n">strain</span> <span class="ow">is</span> <span class="kc">None</span> <span class="ow">and</span> <span class="n">G</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;either strain or G must be present&#39;</span><span class="p">)</span>

    <span class="n">h0</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigHeight</span>
    <span class="n">h</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>

    <span class="nb">print</span> <span class="s1">&#39;Original height </span><span class="si">%.4f</span><span class="s1"> A&#39;</span> <span class="o">%</span> <span class="n">h0</span>
    <span class="nb">print</span> <span class="s1">&#39;Initial measured height </span><span class="si">%.4f</span><span class="s1"> A&#39;</span> <span class="o">%</span> <span class="n">h</span>

    <span class="k">if</span> <span class="n">strain</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">strain</span> <span class="o">=</span> <span class="n">crack_g_to_strain</span><span class="p">(</span><span class="n">G</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="n">eps0</span> <span class="o">=</span> <span class="p">(</span><span class="n">h</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span>
    <span class="n">eps1</span> <span class="o">=</span> <span class="n">strain</span>

    <span class="nb">print</span> <span class="s1">&#39;Initial strain_yy </span><span class="si">%.4f</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">eps0</span>
    <span class="nb">print</span> <span class="s1">&#39;Initial G </span><span class="si">%.4f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span> <span class="n">crack_strain_to_g</span><span class="p">(</span><span class="n">eps0</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">([(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps1</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps0</span><span class="p">),</span> <span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps1</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps0</span><span class="p">),</span> <span class="mf">1.0</span><span class="p">])</span>

    <span class="n">b</span> <span class="o">=</span> <span class="n">transform</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>

    <span class="n">h1</span> <span class="o">=</span> <span class="n">b</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">b</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="nb">print</span> <span class="s1">&#39;Final strain_yy </span><span class="si">% .4f</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="p">((</span><span class="n">h1</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span><span class="p">)</span>

    <span class="n">b</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">b</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">b</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="nb">print</span> <span class="s1">&#39;Final G </span><span class="si">%.4f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span> <span class="n">b</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span>

    <span class="n">crack_update_connect</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">b</span>


<span class="k">def</span> <span class="nf">crack_rescale_uniaxial</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">params</span><span class="p">,</span> <span class="n">strain</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">G</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
    <span class="k">if</span> <span class="n">strain</span> <span class="ow">is</span> <span class="kc">None</span> <span class="ow">and</span> <span class="n">G</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;either strain or G must be present&#39;</span><span class="p">)</span>

    <span class="n">h0</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigHeight</span>
    <span class="n">h</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>

    <span class="nb">print</span> <span class="s1">&#39;Original height </span><span class="si">%.4f</span><span class="s1"> A&#39;</span> <span class="o">%</span> <span class="n">h0</span>
    <span class="nb">print</span> <span class="s1">&#39;Initial measured height </span><span class="si">%.4f</span><span class="s1"> A&#39;</span> <span class="o">%</span> <span class="n">h</span>

    <span class="k">if</span> <span class="n">strain</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">strain</span> <span class="o">=</span> <span class="n">crack_g_to_strain</span><span class="p">(</span><span class="n">G</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="n">eps0</span> <span class="o">=</span> <span class="p">(</span><span class="n">h</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span>
    <span class="n">eps1</span> <span class="o">=</span> <span class="n">strain</span>

    <span class="nb">print</span> <span class="s1">&#39;Initial strain </span><span class="si">%.4f</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">eps0</span>
    <span class="nb">print</span> <span class="s1">&#39;Initial G </span><span class="si">%.4f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span> <span class="n">crack_strain_to_g</span><span class="p">(</span><span class="n">eps0</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">([</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps1</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">+</span><span class="n">eps0</span><span class="p">),</span> <span class="mf">1.0</span><span class="p">])</span>

    <span class="n">b</span> <span class="o">=</span> <span class="n">transform</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">t</span><span class="p">)</span>

    <span class="n">h1</span> <span class="o">=</span> <span class="n">b</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">b</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="nb">print</span> <span class="s1">&#39;Final strain_yy </span><span class="si">% .4f</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="p">((</span><span class="n">h1</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span><span class="p">)</span>

    <span class="n">b</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">b</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">b</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">h0</span><span class="p">)</span>

    <span class="nb">print</span> <span class="s1">&#39;Final G </span><span class="si">%.4f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span> <span class="n">b</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span>

    <span class="n">crack_update_connect</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">params</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">b</span>



<span class="k">def</span> <span class="nf">makecrack_initial_velocity_field</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">,</span> <span class="n">advance_step</span><span class="p">,</span> <span class="n">advance_time</span><span class="p">):</span>
    <span class="n">crack_slab_1</span> <span class="o">=</span> <span class="n">makecrack_main</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">)</span>

    <span class="c1"># advance by one bond</span>
    <span class="n">params</span><span class="o">.</span><span class="n">crack_seed_length</span> <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_seed_length</span> <span class="o">+</span> <span class="n">advance_step</span>
    <span class="n">crack_slab_2</span> <span class="o">=</span> <span class="n">makecrack_main</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">)</span>

    <span class="n">crack_slab_1</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;velo&#39;</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="n">n_cols</span><span class="o">=</span><span class="mi">3</span><span class="p">)</span>
    <span class="n">crack_slab_1</span><span class="o">.</span><span class="n">velo</span><span class="p">[</span><span class="o">...</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">crack_slab_2</span><span class="o">.</span><span class="n">pos</span> <span class="o">-</span> <span class="n">crack_slab_1</span><span class="o">.</span><span class="n">pos</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">advance_step</span><span class="o">*</span><span class="n">advance_time</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">crack_slab_1</span>

<span class="k">def</span> <span class="nf">makecrack</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">):</span>
    <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_initial_velocity_field</span><span class="p">:</span>
        <span class="n">slab</span> <span class="o">=</span> <span class="n">makecrack_initial_velocity_field</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">,</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_initial_velocity_field_dx</span><span class="p">,</span>
                                                <span class="n">params</span><span class="o">.</span><span class="n">crack_initial_velocity_field_dt</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">slab</span> <span class="o">=</span> <span class="n">makecrack_main</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">stem</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">slab</span>


<div class="viewcode-block" id="crack_strain_energy_release_rate"><a class="viewcode-back" href="../../crack.html#quippy.crack.crack_strain_energy_release_rate">[docs]</a><span class="k">def</span> <span class="nf">crack_strain_energy_release_rate</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">bulk</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">f_min</span><span class="o">=.</span><span class="mi">8</span><span class="p">,</span> <span class="n">f_max</span><span class="o">=.</span><span class="mi">9</span><span class="p">,</span> <span class="n">stem</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_pos</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute strain energy release rate G from elastic potential energy in a strip</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="nb">print</span> <span class="s1">&#39;Analytical effective elastic modulus E</span><span class="se">\&#39;</span><span class="s1"> = &#39;</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="o">**</span><span class="mi">2</span><span class="p">),</span> <span class="s1">&#39;GPa&#39;</span>
    <span class="nb">print</span> <span class="s1">&#39;Analytical energy release rate G = &#39;</span><span class="p">,</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigHeight</span><span class="p">),</span> <span class="s1">&#39;J/m^2&#39;</span>

    <span class="k">if</span> <span class="n">bulk</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">stem</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span> <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;Either &quot;bulk&quot; or &quot;stem&quot; must be present&#39;</span><span class="p">)</span>
        <span class="n">bulk</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">stem</span><span class="o">+</span><span class="s1">&#39;_bulk.xyz&#39;</span><span class="p">)</span>

    <span class="k">if</span> <span class="ow">not</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="s1">&#39;local_energy&#39;</span><span class="p">)</span> <span class="ow">or</span> <span class="ow">not</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">bulk</span><span class="p">,</span> <span class="s1">&#39;energy&#39;</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">stem</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span> <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;local_energy property not found in Atoms and &quot;stem&quot; is missing&#39;</span><span class="p">)</span>
        <span class="n">xmlfile</span> <span class="o">=</span> <span class="n">stem</span><span class="o">+</span><span class="s1">&#39;.xml&#39;</span>
        <span class="n">params</span> <span class="o">=</span> <span class="n">CrackParams</span><span class="p">(</span><span class="n">xmlfile</span><span class="p">)</span>
        <span class="n">pot</span> <span class="o">=</span> <span class="n">Potential</span><span class="p">(</span><span class="n">params</span><span class="o">.</span><span class="n">classical_args</span><span class="p">,</span> <span class="n">param_filename</span><span class="o">=</span><span class="n">stem</span><span class="o">+</span><span class="s1">&#39;.xml&#39;</span><span class="p">)</span>
        <span class="n">pot</span><span class="o">.</span><span class="n">print_</span><span class="p">()</span>

        <span class="k">if</span> <span class="ow">not</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="s1">&#39;local_energy&#39;</span><span class="p">):</span>
            <span class="k">if</span> <span class="n">avg_pos</span><span class="p">:</span>
                <span class="n">tmp_pos</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
                <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="o">...</span><span class="p">]</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">avgpos</span>
            <span class="n">at</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">pot</span><span class="o">.</span><span class="n">cutoff</span><span class="p">()</span><span class="o">+</span><span class="mf">1.</span><span class="p">)</span>
            <span class="n">at</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>
            <span class="n">pot</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">at</span><span class="p">,</span> <span class="n">args_str</span><span class="o">=</span><span class="s2">&quot;local_energy&quot;</span><span class="p">)</span>
            <span class="k">if</span> <span class="n">avg_pos</span><span class="p">:</span>
                <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="o">...</span><span class="p">]</span> <span class="o">=</span> <span class="n">tmp_pos</span>

        <span class="k">if</span> <span class="ow">not</span> <span class="nb">hasattr</span><span class="p">(</span><span class="n">bulk</span><span class="p">,</span> <span class="s1">&#39;energy&#39;</span><span class="p">):</span>
            <span class="n">bulk</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">pot</span><span class="o">.</span><span class="n">cutoff</span><span class="p">()</span><span class="o">+</span><span class="mf">1.</span><span class="p">)</span>
            <span class="n">bulk</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>
            <span class="n">pot</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">bulk</span><span class="p">,</span> <span class="n">args_str</span><span class="o">=</span><span class="s1">&#39;energy&#39;</span><span class="p">)</span>

    <span class="n">h</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="n">h0</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigHeight</span>
    <span class="n">strain</span> <span class="o">=</span> <span class="p">(</span><span class="n">h</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span>
    <span class="nb">print</span> <span class="s1">&#39;Applied strain&#39;</span><span class="p">,</span> <span class="n">strain</span>

    <span class="n">x_min</span> <span class="o">=</span> <span class="n">f_min</span><span class="o">*</span><span class="n">at</span><span class="o">.</span><span class="n">OrigWidth</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigWidth</span><span class="o">/</span><span class="mf">2.</span>
    <span class="n">x_max</span> <span class="o">=</span> <span class="n">f_max</span><span class="o">*</span><span class="n">at</span><span class="o">.</span><span class="n">OrigWidth</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigWidth</span><span class="o">/</span><span class="mf">2.</span>
    <span class="n">strip</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">logical_and</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">move_mask</span> <span class="o">==</span> <span class="mi">1</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">logical_and</span><span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">&gt;</span> <span class="n">x_min</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">&lt;</span> <span class="n">x_max</span><span class="p">))</span>
    <span class="n">at</span><span class="o">.</span><span class="n">add_property</span><span class="p">(</span><span class="s1">&#39;strip&#39;</span><span class="p">,</span> <span class="n">strip</span><span class="p">,</span> <span class="n">overwrite</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>

    <span class="n">strip_depth</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">lattice</span><span class="p">[</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">]</span>
    <span class="n">strip_width</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">strip</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">strip</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="n">strip_height</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="n">strip</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="n">strip</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="n">strip_volume</span> <span class="o">=</span> <span class="n">strip_width</span><span class="o">*</span><span class="n">strip_height</span><span class="o">*</span><span class="n">strip_depth</span>
    <span class="nb">print</span> <span class="s1">&#39;Strip contains&#39;</span><span class="p">,</span> <span class="n">strip</span><span class="o">.</span><span class="n">sum</span><span class="p">(),</span> <span class="s1">&#39;atoms&#39;</span><span class="p">,</span> <span class="s1">&#39;width&#39;</span><span class="p">,</span> <span class="n">strip_width</span><span class="p">,</span> <span class="s1">&#39;height&#39;</span><span class="p">,</span> <span class="n">strip_height</span><span class="p">,</span> <span class="s1">&#39;volume&#39;</span><span class="p">,</span> <span class="n">strip_volume</span>

    <span class="n">strain_energy_density</span> <span class="o">=</span> <span class="p">(</span><span class="n">at</span><span class="o">.</span><span class="n">local_energy</span><span class="p">[</span><span class="n">strip</span><span class="p">]</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span> <span class="o">-</span> <span class="n">bulk</span><span class="o">.</span><span class="n">energy</span><span class="o">/</span><span class="n">bulk</span><span class="o">.</span><span class="n">n</span><span class="o">*</span><span class="n">strip</span><span class="o">.</span><span class="n">sum</span><span class="p">())</span><span class="o">/</span><span class="n">strip_volume</span>

    <span class="nb">print</span> <span class="s1">&#39;Strain energy density in strip&#39;</span><span class="p">,</span> <span class="n">strain_energy_density</span><span class="p">,</span> <span class="s1">&#39;eV/A**3&#39;</span>

    <span class="n">E_effective</span> <span class="o">=</span> <span class="mi">2</span><span class="o">*</span><span class="n">strain_energy_density</span><span class="o">/</span><span class="n">strain</span><span class="o">**</span><span class="mi">2</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span>
    <span class="nb">print</span> <span class="s1">&#39;Effective elastic modulus E =&#39;</span><span class="p">,</span> <span class="n">E_effective</span><span class="p">,</span> <span class="s1">&#39;GPa&#39;</span>

    <span class="n">G_effective</span> <span class="o">=</span> <span class="n">strain_energy_density</span><span class="o">*</span><span class="n">strip_height</span><span class="o">*</span><span class="n">J_PER_M2</span>
    <span class="nb">print</span> <span class="s1">&#39;Effective energy release rate G =&#39;</span><span class="p">,</span> <span class="n">G_effective</span><span class="p">,</span> <span class="s1">&#39;J/m^2&#39;</span>

    <span class="k">return</span> <span class="n">G_effective</span></div>


<span class="k">def</span> <span class="nf">crack_make_seed_curved_front</span><span class="p">(</span><span class="n">slab</span><span class="p">,</span> <span class="n">params</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;Given a slab, introduce a 3D curved crack front.&quot;&quot;&quot;</span>

   <span class="n">orig_width</span> <span class="o">=</span> <span class="n">slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigWidth&#39;</span><span class="p">]</span>
   <span class="n">orig_height</span> <span class="o">=</span> <span class="n">slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>
   <span class="n">orig_depth</span>  <span class="o">=</span> <span class="n">slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;OrigDepth&#39;</span><span class="p">]</span>

   <span class="n">crack_length</span> <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_seed_length</span><span class="o">/</span><span class="n">orig_width</span>
   <span class="n">crack_curvature</span> <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_curvature</span>
   <span class="n">crack_depth</span> <span class="o">=</span> <span class="mf">1.0</span>

   <span class="k">if</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_g</span> <span class="o">&gt;</span> <span class="mf">0.</span><span class="p">:</span>
      <span class="n">epsilon_max</span> <span class="o">=</span> <span class="n">crack_g_to_strain</span><span class="p">(</span><span class="n">params</span><span class="o">.</span><span class="n">crack_g</span><span class="p">,</span>
                                     <span class="n">slab</span><span class="o">.</span><span class="n">YoungsModulus</span><span class="p">,</span>
                                     <span class="n">slab</span><span class="o">.</span><span class="n">PoissonRatio_yx</span><span class="p">,</span>
                                     <span class="n">orig_height</span><span class="p">)</span>
   <span class="k">else</span><span class="p">:</span>
      <span class="n">epsilon_max</span> <span class="o">=</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_strain</span>

   <span class="c1"># ... then shift so coordinates are in ranges x=[0..a], y=[0..b], z=[0..c]</span>
   <span class="n">slab</span><span class="o">.</span><span class="n">pos</span> <span class="o">+=</span> <span class="n">np</span><span class="o">.</span><span class="n">tile</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">orig_width</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">,</span> <span class="n">orig_depth</span><span class="p">])[:,</span><span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">]</span><span class="o">/</span><span class="mf">2.</span><span class="p">,</span> <span class="n">slab</span><span class="o">.</span><span class="n">n</span><span class="p">)</span>

   <span class="n">xlin</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">crack_length</span><span class="o">*</span><span class="n">orig_width</span>
   <span class="n">crack_center</span> <span class="o">=</span> <span class="n">farray</span><span class="p">([</span><span class="n">crack_length</span><span class="o">*</span><span class="n">orig_width</span><span class="p">,</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">orig_height</span><span class="p">,</span> <span class="mf">0.</span><span class="p">])</span>

   <span class="n">epsilon_min</span> <span class="o">=</span> <span class="mf">0.</span>
   <span class="n">epsilon_local</span> <span class="o">=</span> <span class="mf">0.</span>
   <span class="n">dy</span> <span class="o">=</span> <span class="mf">0.</span>

   <span class="n">a</span> <span class="o">=</span> <span class="n">crack_curvature</span>
   <span class="n">x1</span> <span class="o">=</span> <span class="n">crack_length</span><span class="o">*</span><span class="n">orig_width</span>
   <span class="n">x2</span> <span class="o">=</span> <span class="n">x1</span> <span class="o">+</span> <span class="n">crack_depth</span>
   <span class="n">z1</span> <span class="o">=</span> <span class="mf">0.</span>
   <span class="n">z2</span> <span class="o">=</span> <span class="n">orig_depth</span>
   <span class="n">slope</span> <span class="o">=</span> <span class="p">(</span><span class="n">x2</span><span class="o">-</span><span class="n">x1</span><span class="p">)</span><span class="o">/</span><span class="p">(</span><span class="n">z2</span><span class="o">-</span><span class="n">z1</span><span class="p">)</span>

   <span class="n">dymax</span> <span class="o">=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">epsilon_max</span><span class="o">*</span><span class="n">orig_height</span>

   <span class="n">z</span> <span class="o">=</span> <span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">3</span><span class="p">,:]</span>
   <span class="n">crack_x_center</span> <span class="o">=</span> <span class="n">crack_curvature</span><span class="o">*</span><span class="p">(</span><span class="n">z</span><span class="o">*</span><span class="n">z</span> <span class="o">-</span> <span class="n">z1</span><span class="o">*</span><span class="n">z1</span><span class="p">)</span> <span class="o">+</span> <span class="p">(</span><span class="n">slope</span><span class="o">-</span><span class="n">crack_curvature</span><span class="o">*</span><span class="p">(</span><span class="n">z2</span><span class="o">+</span><span class="n">z1</span><span class="p">))</span><span class="o">*</span><span class="p">(</span><span class="n">z</span><span class="o">-</span><span class="n">z1</span><span class="p">)</span> <span class="o">+</span> <span class="n">crack_center</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span>
   <span class="c1">##slab.add_property(&#39;xc&#39;, crack_x_center, overwrite=True)</span>

   <span class="n">dymin</span> <span class="o">=</span> <span class="nb">abs</span><span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">-</span><span class="n">crack_center</span><span class="p">[</span><span class="mi">2</span><span class="p">])</span><span class="o">*</span><span class="n">epsilon_max</span>

   <span class="n">dy</span> <span class="o">=</span> <span class="n">dymin</span>
   <span class="n">dy</span><span class="p">[</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">&lt;</span> <span class="n">xlin</span><span class="p">]</span> <span class="o">=</span> <span class="n">dymax</span>
   <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">&gt;</span> <span class="n">xlin</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,:]</span> <span class="o">&lt;=</span> <span class="n">crack_x_center</span><span class="p">)</span>
   <span class="c1">##slab.add_property(&#39;mask&#39;, mask, overwrite=True)</span>
   <span class="n">dy</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">dymax</span> <span class="o">-</span> <span class="n">dymin</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span><span class="o">/</span><span class="p">(</span><span class="n">xlin</span> <span class="o">-</span> <span class="n">crack_x_center</span><span class="p">[</span><span class="n">mask</span><span class="p">])</span><span class="o">*</span><span class="nb">abs</span><span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">1</span><span class="p">,</span><span class="n">mask</span><span class="p">]</span> <span class="o">-</span> <span class="n">xlin</span><span class="p">)</span> <span class="o">+</span> <span class="n">dymax</span>

   <span class="n">dy</span> <span class="o">*=</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span> <span class="o">-</span> <span class="n">crack_center</span><span class="p">[</span><span class="mi">2</span><span class="p">])</span>
   <span class="c1">##slab.add_property(&#39;dy&#39;, dy, overwrite=True)</span>
   <span class="n">slab</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span> <span class="o">+=</span> <span class="n">dy</span>

   <span class="c1"># Centre slab in cell</span>
   <span class="n">slab</span><span class="o">.</span><span class="n">pos</span> <span class="o">+=</span> <span class="n">np</span><span class="o">.</span><span class="n">tile</span><span class="p">((</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span><span class="o">/</span><span class="mf">2.</span><span class="p">]</span><span class="o">*</span><span class="mi">3</span><span class="p">)</span> <span class="o">+</span>
                        <span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">slab</span><span class="o">.</span><span class="n">lattice</span><span class="p">)</span><span class="o">/</span><span class="mf">2.</span><span class="p">)[:,</span><span class="n">np</span><span class="o">.</span><span class="n">newaxis</span><span class="p">],</span> <span class="n">slab</span><span class="o">.</span><span class="n">n</span><span class="p">)</span>
   <span class="n">slab</span><span class="o">.</span><span class="n">map_into_cell</span><span class="p">()</span>

   <span class="n">slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;CrackPosx&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_center</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span><span class="o">/</span><span class="mf">2.</span>
   <span class="n">slab</span><span class="o">.</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;CrackPosy&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">crack_center</span><span class="p">[</span><span class="mi">2</span><span class="p">]</span> <span class="o">+</span> <span class="n">params</span><span class="o">.</span><span class="n">crack_vacuum_size</span><span class="o">/</span><span class="mf">2.</span>


<div class="viewcode-block" id="crack_strain"><a class="viewcode-back" href="../../crack.html#quippy.crack.crack_strain">[docs]</a><span class="k">def</span> <span class="nf">crack_strain</span><span class="p">(</span><span class="n">at</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   Returns strain of crack specimen</span>
<span class="sd">   &quot;&quot;&quot;</span>
   <span class="n">h</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="mi">2</span><span class="p">,:]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
   <span class="n">h0</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">OrigHeight</span>
   <span class="n">eps</span> <span class="o">=</span> <span class="p">(</span><span class="n">h</span> <span class="o">-</span> <span class="n">h0</span><span class="p">)</span><span class="o">/</span><span class="n">h0</span>
   <span class="k">return</span> <span class="n">eps</span></div>

<span class="k">def</span> <span class="nf">energy_difference</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">b</span><span class="p">,</span><span class="n">pot</span><span class="p">,</span><span class="n">eps</span><span class="p">,</span><span class="n">relax</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   Compute energy difference between configurations `a` and `b` using</span>
<span class="sd">   `pot`, after the application of a homogenous-xy rescaling to strain</span>
<span class="sd">   `eps`. If `relax`=True, the rescaled confiugrations are</span>
<span class="sd">   structurally optimised.</span>
<span class="sd">   &quot;&quot;&quot;</span>

   <span class="n">ap</span> <span class="o">=</span> <span class="n">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">cp</span><span class="p">,</span> <span class="n">eps</span><span class="p">)</span>
   <span class="n">bp</span> <span class="o">=</span> <span class="n">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">cp</span><span class="p">,</span> <span class="n">eps</span><span class="p">)</span>
   <span class="n">ap</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">pot</span><span class="o">.</span><span class="n">cutoff</span><span class="p">()</span><span class="o">+</span><span class="mf">1.0</span><span class="p">)</span>
   <span class="n">ap</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>
   <span class="n">bp</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">pot</span><span class="o">.</span><span class="n">cutoff</span><span class="p">()</span><span class="o">+</span><span class="mf">1.0</span><span class="p">)</span>
   <span class="n">bp</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>
   <span class="k">if</span> <span class="n">relax</span><span class="p">:</span>
      <span class="n">pot</span><span class="o">.</span><span class="n">minim</span><span class="p">(</span><span class="n">ap</span><span class="p">,</span> <span class="s1">&#39;cg&#39;</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mi">100</span><span class="p">,</span> <span class="n">do_pos</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">do_lat</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
      <span class="n">pot</span><span class="o">.</span><span class="n">minim</span><span class="p">(</span><span class="n">bp</span><span class="p">,</span> <span class="s1">&#39;cg&#39;</span><span class="p">,</span> <span class="mf">0.1</span><span class="p">,</span> <span class="mi">100</span><span class="p">,</span> <span class="n">do_pos</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span> <span class="n">do_lat</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
   <span class="n">pot</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">ap</span><span class="p">,</span> <span class="n">args_str</span><span class="o">=</span><span class="s2">&quot;energy&quot;</span><span class="p">)</span>
   <span class="n">pot</span><span class="o">.</span><span class="n">calc</span><span class="p">(</span><span class="n">bp</span><span class="p">,</span> <span class="n">args_str</span><span class="o">=</span><span class="s2">&quot;energy&quot;</span><span class="p">)</span>
   <span class="nb">print</span> <span class="s2">&quot;ENERGY_DIFFERENCE </span><span class="si">%f</span><span class="s2"> </span><span class="si">%f</span><span class="s2"> </span><span class="si">%f</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">crack_strain</span><span class="p">(</span><span class="n">ap</span><span class="p">),</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">ap</span><span class="p">),</span> <span class="n">ap</span><span class="o">.</span><span class="n">energy</span><span class="o">-</span><span class="n">bp</span><span class="o">.</span><span class="n">energy</span><span class="p">)</span>
   <span class="k">return</span> <span class="n">ap</span><span class="o">.</span><span class="n">energy</span><span class="o">-</span><span class="n">bp</span><span class="o">.</span><span class="n">energy</span>


<div class="viewcode-block" id="crack_find_griffith_load"><a class="viewcode-back" href="../../crack.html#quippy.crack.crack_find_griffith_load">[docs]</a><span class="k">def</span> <span class="nf">crack_find_griffith_load</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">,</span> <span class="n">pot</span><span class="p">,</span> <span class="n">relax</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   Given two configurations (a, b) which differ by one broken bond,</span>
<span class="sd">   find the Griffith load, that is the load at which `a` and `b` have</span>
<span class="sd">   the same energy accorinding to the model potential `pot`.</span>

<span class="sd">   Returns (strain, G, a_rescaled, b_rescaled).</span>
<span class="sd">   &quot;&quot;&quot;</span>

   <span class="n">eps0</span> <span class="o">=</span> <span class="n">crack_strain</span><span class="p">(</span><span class="n">a</span><span class="p">)</span>

   <span class="k">def</span> <span class="nf">func</span><span class="p">(</span><span class="n">x</span><span class="p">):</span>
       <span class="n">eps</span><span class="p">,</span> <span class="o">=</span> <span class="n">x</span>
       <span class="k">return</span> <span class="n">energy_difference</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">b</span><span class="p">,</span><span class="n">p</span><span class="p">,</span><span class="n">eps</span><span class="p">,</span><span class="n">relax</span><span class="o">=</span><span class="n">relax</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span>

   <span class="n">ds</span> <span class="o">=</span> <span class="n">DownhillSimplex</span><span class="p">(</span><span class="n">func</span><span class="p">,</span> <span class="n">x0</span><span class="o">=</span><span class="p">[</span><span class="n">eps0</span><span class="p">],</span> <span class="n">deltas</span><span class="o">=</span><span class="p">[</span><span class="mf">0.001</span><span class="p">],</span> <span class="n">ftol</span><span class="o">=</span><span class="mf">0.01</span><span class="p">)</span>
   <span class="n">eps</span><span class="p">,</span> <span class="o">=</span> <span class="n">ds</span><span class="o">.</span><span class="n">minimise</span><span class="p">()</span>
   <span class="n">ap</span> <span class="o">=</span> <span class="n">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">cp</span><span class="p">,</span> <span class="n">eps</span><span class="p">)</span>
   <span class="n">bp</span> <span class="o">=</span> <span class="n">crack_rescale_homogeneous_xy</span><span class="p">(</span><span class="n">b</span><span class="p">,</span> <span class="n">cp</span><span class="p">,</span> <span class="n">eps</span><span class="p">)</span>

   <span class="nb">print</span> <span class="s1">&#39;Griffith critical strain = </span><span class="si">%f</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">crack_strain</span><span class="p">(</span><span class="n">ap</span><span class="p">)</span>
   <span class="nb">print</span> <span class="s1">&#39;Griffith critical G = </span><span class="si">%f</span><span class="s1"> J/m^2&#39;</span> <span class="o">%</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">ap</span><span class="p">)</span>

   <span class="k">return</span> <span class="p">(</span><span class="n">crack_strain</span><span class="p">(</span><span class="n">ap</span><span class="p">),</span> <span class="n">crack_measure_g</span><span class="p">(</span><span class="n">ap</span><span class="p">),</span> <span class="n">ap</span><span class="p">,</span> <span class="n">bp</span><span class="p">)</span></div>


<span class="k">def</span> <span class="nf">find_mapping</span><span class="p">(</span><span class="n">at_in</span><span class="p">,</span> <span class="n">vectors</span><span class="p">,</span> <span class="n">lattice_constant</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-4</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   Find a mapping between pairs of atoms displaced by a given vector</span>
<span class="sd">   (or by one of a number of given vectors).</span>
<span class="sd">   &quot;&quot;&quot;</span>

   <span class="k">class</span> <span class="nc">FoundMapping</span><span class="p">:</span>
      <span class="k">pass</span>

   <span class="n">mapping</span> <span class="o">=</span> <span class="n">fzeros</span><span class="p">(</span><span class="n">at_in</span><span class="o">.</span><span class="n">n</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">int32</span><span class="p">)</span>
   <span class="n">at</span> <span class="o">=</span> <span class="n">at_in</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

   <span class="c1"># cutoff should be larger than all displacment vectors</span>
   <span class="n">vlen</span> <span class="o">=</span> <span class="p">[</span><span class="n">farray</span><span class="p">(</span><span class="n">vector</span><span class="p">)</span><span class="o">.</span><span class="n">norm</span><span class="p">()</span><span class="o">*</span><span class="n">lattice_constant</span> <span class="k">for</span> <span class="n">vector</span> <span class="ow">in</span> <span class="n">vectors</span><span class="p">]</span>
   <span class="n">at</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="nb">max</span><span class="p">(</span><span class="n">vlen</span><span class="p">)</span><span class="o">+</span><span class="mf">0.5</span><span class="p">)</span>
   <span class="nb">print</span> <span class="s1">&#39;cutoff = &#39;</span><span class="p">,</span> <span class="n">at</span><span class="o">.</span><span class="n">cutoff</span>
   <span class="n">at</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>

   <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">indices</span><span class="p">:</span>
      <span class="k">try</span><span class="p">:</span>
         <span class="k">for</span> <span class="n">neighb</span> <span class="ow">in</span> <span class="n">at</span><span class="o">.</span><span class="n">neighbours</span><span class="p">[</span><span class="n">i</span><span class="p">]:</span>
            <span class="k">for</span> <span class="n">vector</span> <span class="ow">in</span> <span class="n">vectors</span><span class="p">:</span>
               <span class="k">if</span> <span class="p">((</span><span class="n">neighb</span><span class="o">.</span><span class="n">diff</span><span class="o">/</span><span class="n">lattice_constant</span> <span class="o">-</span> <span class="n">vector</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">:</span>
                  <span class="nb">print</span> <span class="n">i</span><span class="p">,</span> <span class="n">neighb</span><span class="o">.</span><span class="n">j</span><span class="p">,</span> <span class="n">vector</span>
                  <span class="n">mapping</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">neighb</span><span class="o">.</span><span class="n">j</span>
                  <span class="k">raise</span> <span class="n">FoundMapping</span>

         <span class="c1"># corresponding atom is off the edge, so map atom onto itself</span>
         <span class="nb">print</span> <span class="n">i</span><span class="p">,</span> <span class="s1">&#39;not found!&#39;</span>
         <span class="n">mapping</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">i</span>
      <span class="k">except</span> <span class="n">FoundMapping</span><span class="p">:</span>
         <span class="k">continue</span>

   <span class="k">return</span> <span class="n">mapping</span>


<span class="n">crack_advance_table</span> <span class="o">=</span> <span class="p">{</span>
   <span class="s1">&#39;Si(110)[001b]+1&#39;</span><span class="p">:</span>  <span class="p">([[</span><span class="o">-</span><span class="mf">1.0</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">)),</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.25</span><span class="p">],</span>
                        <span class="p">[</span><span class="o">-</span><span class="mf">1.0</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">)),</span> <span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">0.25</span><span class="p">]],</span>
                        <span class="p">[[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">1.0</span><span class="p">]]),</span>
   <span class="s1">&#39;Si(110)[001b]-1&#39;</span><span class="p">:</span>  <span class="p">([[</span><span class="mf">1.0</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">)),</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.25</span><span class="p">],</span>
                        <span class="p">[</span><span class="mf">1.0</span><span class="o">/</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">)),</span> <span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">0.25</span><span class="p">]],</span>
                        <span class="p">[[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">1.0</span><span class="p">]]),</span>
   <span class="s1">&#39;Si(110)[001b]+2&#39;</span><span class="p">:</span>  <span class="p">([[</span><span class="o">-</span><span class="mf">1.0</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">),</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">]],</span>
                        <span class="p">[[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">]]),</span>
   <span class="s1">&#39;Si(110)[001b]-2&#39;</span><span class="p">:</span>  <span class="p">([[</span><span class="mf">1.0</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="p">),</span> <span class="mf">0.0</span><span class="p">,</span>  <span class="mf">0.0</span><span class="p">]],</span>
                        <span class="p">[[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">],</span>
                         <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">]])</span>
<span class="p">}</span>


<span class="k">def</span> <span class="nf">crack_unit_advance</span><span class="p">(</span><span class="n">slab</span><span class="p">,</span> <span class="n">at</span><span class="p">,</span> <span class="n">lattice_constant</span><span class="p">,</span> <span class="n">vectors</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">transformation</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">system</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
   <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">   Given isoatomic slab and crack configurations, return copy with crack advanced by one bond</span>
<span class="sd">   &quot;&quot;&quot;</span>

   <span class="k">if</span> <span class="n">system</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
      <span class="k">try</span><span class="p">:</span>
         <span class="n">vectors</span><span class="p">,</span> <span class="n">transformation</span> <span class="o">=</span> <span class="n">crack_advance_table</span><span class="p">[</span><span class="n">system</span><span class="p">]</span>
      <span class="k">except</span> <span class="ne">KeyError</span><span class="p">:</span>
         <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;unknown crack system </span><span class="si">%s</span><span class="s2">&quot;</span> <span class="o">%</span> <span class="n">system</span><span class="p">)</span>

   <span class="k">if</span> <span class="n">vectors</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
      <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s2">&quot;missing vectors&quot;</span><span class="p">)</span>

   <span class="c1"># find mapping from old to new atom indices</span>
   <span class="n">mapping</span> <span class="o">=</span> <span class="n">find_mapping</span><span class="p">(</span><span class="n">slab</span><span class="p">,</span> <span class="n">vectors</span><span class="p">,</span> <span class="n">lattice_constant</span><span class="p">)</span>

   <span class="n">disp</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">pos</span> <span class="o">-</span> <span class="n">slab</span><span class="o">.</span><span class="n">pos</span>
   <span class="n">new_disp</span> <span class="o">=</span> <span class="n">disp</span><span class="p">[</span><span class="n">mapping</span><span class="p">]</span>

   <span class="c1"># optionally apply a transformation to the new displacements</span>
   <span class="k">if</span> <span class="n">transformation</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
      <span class="n">new_disp</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">transformation</span><span class="p">,</span> <span class="n">new_disp</span><span class="p">)</span>

   <span class="n">shift</span> <span class="o">=</span> <span class="p">[</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">,</span> <span class="n">lattice_constant</span><span class="o">/</span><span class="mf">2.0</span><span class="p">]</span>

   <span class="n">result</span> <span class="o">=</span> <span class="n">at</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
   <span class="n">result</span><span class="o">.</span><span class="n">pos</span><span class="p">[</span><span class="o">...</span><span class="p">]</span> <span class="o">=</span> <span class="n">slab</span><span class="o">.</span><span class="n">pos</span> <span class="o">+</span> <span class="n">new_disp</span>
   <span class="k">return</span> <span class="n">result</span>


<span class="c1"># displacement vectors from an atom to the atom horizontally ahead</span>
<span class="c1"># of it (along x), for a few different crack systems, in scaled</span>
<span class="c1"># coordinates. There are two possibilities in each case, one in</span>
<span class="c1"># front and one behind (along z axis)</span>

<span class="n">crack_advance_step_frac_coords</span> <span class="o">=</span> <span class="p">{</span>
    <span class="s1">&#39;(111)[01-1]&#39;</span><span class="p">:</span>
        <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">6.</span><span class="p">)</span> <span class="o">/</span> <span class="mi">4</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">+</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mi">4</span><span class="p">],</span>
                  <span class="p">[</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">6.</span><span class="p">)</span> <span class="o">/</span> <span class="mi">4</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">-</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mf">4.</span><span class="p">]]),</span>

    <span class="s1">&#39;(110)[001]&#39;</span><span class="p">:</span>
        <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mf">4.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="mf">0.25</span><span class="p">],</span>
                  <span class="p">[</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mf">4.</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">-</span><span class="mf">0.25</span><span class="p">]]),</span>

    <span class="s1">&#39;(110)[1-10]&#39;</span><span class="p">:</span>
        <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([[</span><span class="mf">0.25</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">+</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mi">4</span><span class="p">],</span>
                  <span class="p">[</span><span class="mf">0.75</span><span class="p">,</span> <span class="mf">0.</span><span class="p">,</span> <span class="o">-</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.</span><span class="p">)</span> <span class="o">/</span> <span class="mi">4</span><span class="p">]])</span>
    <span class="p">}</span>


<div class="viewcode-block" id="make_crack_advance_map"><a class="viewcode-back" href="../../crack.html#quippy.crack.make_crack_advance_map">[docs]</a><span class="k">def</span> <span class="nf">make_crack_advance_map</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">tol</span><span class="o">=</span><span class="mf">1e-3</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Find mapping from atom indices to the index of atom one step ahead</span>
<span class="sd">    of them in the crack propagation direction (i.e. along +x).</span>

<span class="sd">    Requires &#39;LatticeConstant&#39;, &#39;CleavagePlane&#39;, and &#39;CrackFront&#39; to</span>
<span class="sd">    be available in atoms.info dictionary.</span>

<span class="sd">    Returns integer array of shape (len(atoms),), and also adds a new</span>
<span class="sd">    array &#39;advance_map&#39; into the Atoms object.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">find_atom_with_displacement</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">target_disps</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Return index of atom with relative displacement from i in target_disps</span>

<span class="sd">        Requires connnectivity to be calculated already</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">indices</span><span class="p">,</span> <span class="n">offsets</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">neighbours</span><span class="o">.</span><span class="n">get_neighbors</span><span class="p">(</span><span class="n">i</span><span class="p">)</span>
        <span class="n">diffs</span> <span class="o">=</span> <span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[</span><span class="n">indices</span><span class="p">]</span> <span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">dot</span><span class="p">(</span><span class="n">offsets</span><span class="p">,</span> <span class="n">atoms</span><span class="o">.</span><span class="n">cell</span><span class="p">)</span> <span class="o">-</span>
                 <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
        <span class="k">for</span> <span class="n">j</span><span class="p">,</span> <span class="n">diff</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">indices</span><span class="p">,</span> <span class="n">diffs</span><span class="p">):</span>
            <span class="k">for</span> <span class="n">target</span> <span class="ow">in</span> <span class="n">target_disps</span><span class="p">:</span>
                <span class="k">if</span> <span class="nb">all</span><span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">diff</span> <span class="o">-</span> <span class="n">target</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">tol</span><span class="p">):</span>
                    <span class="k">return</span> <span class="n">j</span>
        <span class="k">return</span> <span class="mi">0</span>

    <span class="n">a0</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;LatticeConstant&#39;</span><span class="p">]</span>
    <span class="n">cleavage_plane</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CleavagePlane&#39;</span><span class="p">]</span>
    <span class="n">crack_front</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackFront&#39;</span><span class="p">]</span>

    <span class="c1"># lookup this crack system in the dictionary of crack advance steps</span>
    <span class="n">cleavage_plane</span> <span class="o">=</span> <span class="n">MillerPlane</span><span class="p">(</span><span class="n">cleavage_plane</span><span class="p">)</span>
    <span class="n">crack_front</span> <span class="o">=</span> <span class="n">MillerDirection</span><span class="p">(</span><span class="n">crack_front</span><span class="p">)</span>
    <span class="n">key</span> <span class="o">=</span> <span class="nb">str</span><span class="p">(</span><span class="n">cleavage_plane</span><span class="p">)</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="n">crack_front</span><span class="p">)</span>
    <span class="n">steps</span> <span class="o">=</span> <span class="n">a0</span> <span class="o">*</span> <span class="n">crack_advance_step_frac_coords</span><span class="p">[</span><span class="n">key</span><span class="p">]</span>

    <span class="n">advance_map</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">),</span> <span class="n">dtype</span><span class="o">=</span><span class="nb">int</span><span class="p">)</span>

    <span class="c1"># find biggest distance between atoms before and after a crack advance step</span>
    <span class="n">max_step_length</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">steps</span> <span class="o">**</span> <span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">))</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

    <span class="c1"># convert from ase.Atoms to quippy.Atoms, so we can use faster</span>
    <span class="c1"># neighbour lists in Fortran code</span>
    <span class="n">tmp_atoms</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">set_cutoff</span><span class="p">(</span><span class="n">max_step_length</span> <span class="o">+</span> <span class="o">.</span><span class="mi">1</span><span class="p">)</span>
    <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>

    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">tmp_atoms</span><span class="p">)):</span>
        <span class="n">advance_map</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">find_atom_with_displacement</span><span class="p">(</span><span class="n">tmp_atoms</span><span class="p">,</span> <span class="n">i</span><span class="p">,</span> <span class="n">steps</span><span class="p">)</span>

    <span class="c1"># save the map inside the Atoms object, and return a copy</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">new_array</span><span class="p">(</span><span class="s1">&#39;advance_map&#39;</span><span class="p">,</span> <span class="n">advance_map</span><span class="p">)</span>
    <span class="k">return</span> <span class="n">advance_map</span></div>


<div class="viewcode-block" id="find_crack_tip_coordination"><a class="viewcode-back" href="../../crack.html#quippy.crack.find_crack_tip_coordination">[docs]</a><span class="k">def</span> <span class="nf">find_crack_tip_coordination</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">edge_tol</span><span class="o">=</span><span class="mf">10.0</span><span class="p">,</span>
                                <span class="n">strip_height</span><span class="o">=</span><span class="mf">30.0</span><span class="p">,</span> <span class="n">nneightol</span><span class="o">=</span><span class="mf">1.3</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return position of crack tip in `atoms`, based on atomic coordination.</span>

<span class="sd">    If `atoms` does not contain an `advance_map` property, then</span>
<span class="sd">    :func:`make_crack_advance_map` is called to generate the map.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    atoms : :class:`~.Atoms&#39; object</span>
<span class="sd">       The Atoms object containing the crack slab.</span>
<span class="sd">    edge_tol : float</span>
<span class="sd">       Distance from edge of system within which to exclude</span>
<span class="sd">       undercoodinated atoms.</span>
<span class="sd">    strip_height : float</span>
<span class="sd">       Height of strip along centre of slab in which to look</span>
<span class="sd">       for the track.</span>
<span class="sd">    nneightol : float</span>
<span class="sd">       Nearest neighbour tolerance, as a fraction of sum of</span>
<span class="sd">       covalent radii of atomic species.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    crack_pos : array</span>
<span class="sd">       x, y, and z coordinates of the crack tip. Also set in ``CrackPos``</span>
<span class="sd">       in ``atoms.info`` dictionary.</span>
<span class="sd">    tip_atoms : array</span>
<span class="sd">       Indices of atoms near the tip Also set in ``crack_tip`` property.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">old_tip_pos_y</span> <span class="o">=</span> <span class="mi">0</span>
    <span class="k">if</span> <span class="s1">&#39;CrackPos&#39;</span> <span class="ow">in</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">:</span>
        <span class="n">old_tip_pos_y</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span>

    <span class="c1"># Make a copy of atoms as a quippy.Atoms instance, overwriting</span>
    <span class="c1"># positions with time-averages values if they are available, and</span>
    <span class="c1"># then calculate connectivity using nneightol</span>
    <span class="n">tmp_atoms</span> <span class="o">=</span> <span class="n">Atoms</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="k">if</span> <span class="s1">&#39;avgpos&#39;</span> <span class="ow">in</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">arrays</span><span class="p">:</span>
        <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">set_positions</span><span class="p">(</span><span class="n">tmp_atoms</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;avgpos&#39;</span><span class="p">])</span>
    <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">calc_connect</span><span class="p">()</span>

    <span class="n">nn</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">n_neighb</span>
    <span class="n">x</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span>

    <span class="c1"># find undercoordinated atoms in a central strip, and not too</span>
    <span class="c1"># close to the left or right edges</span>
    <span class="n">left</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="n">right</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>
    <span class="n">uc</span> <span class="o">=</span> <span class="p">((</span><span class="n">nn</span> <span class="o">&lt;</span> <span class="mi">4</span><span class="p">)</span> <span class="o">&amp;</span>
          <span class="p">(</span><span class="nb">abs</span><span class="p">(</span><span class="n">y</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">strip_height</span><span class="p">)</span> <span class="o">&amp;</span>
          <span class="p">(</span><span class="n">x</span> <span class="o">&gt;</span> <span class="n">left</span> <span class="o">+</span> <span class="n">edge_tol</span><span class="p">)</span> <span class="o">&amp;</span>
          <span class="p">(</span><span class="n">x</span> <span class="o">&lt;</span> <span class="n">right</span> <span class="o">-</span> <span class="n">edge_tol</span><span class="p">))</span>

    <span class="c1"># position of furthest forward undercoordinated atom ABOVE old tip position</span>
    <span class="n">x_above</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">uc</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">y</span> <span class="o">&gt;</span> <span class="n">old_tip_pos_y</span><span class="p">)]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

    <span class="c1"># position of furthest forward undercoordinated atom BELOW old tip position</span>
    <span class="n">x_below</span> <span class="o">=</span> <span class="n">x</span><span class="p">[</span><span class="n">uc</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">y</span> <span class="o">&lt;</span> <span class="n">old_tip_pos_y</span><span class="p">)]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span>

    <span class="c1"># rightmost undercoordinated atoms, both above and below old tip</span>
    <span class="n">rightmost_uc</span> <span class="o">=</span> <span class="n">uc</span> <span class="o">&amp;</span> <span class="p">(((</span><span class="n">y</span> <span class="o">&gt;</span> <span class="n">old_tip_pos_y</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">x_above</span> <span class="o">==</span> <span class="n">x</span><span class="p">))</span> <span class="o">|</span>
                         <span class="p">((</span><span class="n">y</span> <span class="o">&lt;=</span> <span class="n">old_tip_pos_y</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">x_below</span> <span class="o">==</span> <span class="n">x</span><span class="p">)))</span>

    <span class="c1"># we want the NEXT pair of atoms, so we use the saved mapping from</span>
    <span class="c1"># atom indices to the indices of atoms one unit cell to the right</span>
    <span class="k">if</span> <span class="s1">&#39;advance_map&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">atoms</span><span class="o">.</span><span class="n">arrays</span><span class="p">:</span>
        <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Generating crack advance map...&#39;</span><span class="p">)</span>
        <span class="n">make_crack_advance_map</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>

    <span class="n">advance_map</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;advance_map&#39;</span><span class="p">]</span>
    <span class="n">tip_atoms</span> <span class="o">=</span> <span class="n">advance_map</span><span class="p">[</span><span class="n">rightmost_uc</span><span class="p">]</span>
    <span class="n">tip_pos</span> <span class="o">=</span> <span class="n">tmp_atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[</span><span class="n">tip_atoms</span><span class="p">,</span> <span class="p">:]</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">)</span>

    <span class="c1"># Also save results in Atoms (useful for visualisation)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">tip_pos</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">set_array</span><span class="p">(</span><span class="s1">&#39;crack_tip&#39;</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="kc">False</span><span class="p">]</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">atoms</span><span class="p">)))</span>
    <span class="n">crack_tip</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">arrays</span><span class="p">[</span><span class="s1">&#39;crack_tip&#39;</span><span class="p">]</span>
    <span class="n">crack_tip</span><span class="p">[</span><span class="n">tip_atoms</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
    <span class="k">return</span> <span class="n">tip_pos</span></div>


<div class="viewcode-block" id="irwin_modeI_crack_tip_stress_field"><a class="viewcode-back" href="../../crack.html#quippy.crack.irwin_modeI_crack_tip_stress_field">[docs]</a><span class="k">def</span> <span class="nf">irwin_modeI_crack_tip_stress_field</span><span class="p">(</span><span class="n">K</span><span class="p">,</span> <span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="n">xy_only</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span>
                                       <span class="n">nu</span><span class="o">=</span><span class="mf">0.5</span><span class="p">,</span> <span class="n">stress_state</span><span class="o">=</span><span class="s1">&#39;plane strain&#39;</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Compute Irwin singular crack tip stress field</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    K : float</span>
<span class="sd">       Mode I stress intensity factor. Units should match units of `r`.</span>
<span class="sd">    r : array_like</span>
<span class="sd">       Radial distances from crack tip. Can be a multidimensional</span>
<span class="sd">       array to evaluate stress field on a grid.</span>
<span class="sd">    t : array_like</span>
<span class="sd">       Angles from horzontal line y=0 ahead of crack tip,</span>
<span class="sd">       measured anticlockwise. Should have same shape as `r`.</span>
<span class="sd">    xy_only : bool</span>
<span class="sd">       If True (default) only xx, yy, xy and yx components will be set.</span>
<span class="sd">    nu : float</span>
<span class="sd">       Poisson ratio. Used only when ``xy_only=False``, to determine zz stresses</span>
<span class="sd">    stress_state : str</span>
<span class="sd">       One of&quot;plane stress&quot; or &quot;plane strain&quot;. Used if xyz_only=False to</span>
<span class="sd">       determine zz stresses.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    sigma : array with shape ``r.shape + (3,3)``</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">if</span> <span class="n">r</span><span class="o">.</span><span class="n">shape</span> <span class="o">!=</span> <span class="n">t</span><span class="o">.</span><span class="n">shape</span><span class="p">:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;shapes of radial and angular arrays &quot;r&quot; and &quot;t&quot; must match&#39;</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">stress_state</span> <span class="ow">not</span> <span class="ow">in</span> <span class="p">[</span><span class="s1">&#39;plane strain&#39;</span><span class="p">,</span> <span class="s1">&#39;plane stress&#39;</span><span class="p">]:</span>
        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;stress_state should be either &quot;plane strain&quot; or &quot;plane stress&quot;.&#39;</span><span class="p">)</span>

    <span class="n">sigma</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span><span class="n">r</span><span class="o">.</span><span class="n">shape</span> <span class="o">+</span> <span class="p">(</span><span class="mi">3</span><span class="p">,</span> <span class="mi">3</span><span class="p">))</span>
    <span class="n">radial</span> <span class="o">=</span> <span class="n">K</span><span class="o">*</span><span class="mf">1.</span><span class="o">/</span><span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">pi</span><span class="o">*</span><span class="n">r</span><span class="p">)</span>

    <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">radial</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="mf">1.0</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">))</span> <span class="c1"># xx</span>
    <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">radial</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="p">(</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">))</span> <span class="c1"># yy</span>
    <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">radial</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">*</span><span class="n">np</span><span class="o">.</span><span class="n">cos</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">t</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span>         <span class="c1"># xy</span>
    <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span>                                               <span class="c1"># yx=xy</span>

    <span class="k">if</span> <span class="ow">not</span> <span class="n">xy_only</span> <span class="ow">and</span> <span class="n">stress_state</span> <span class="o">==</span> <span class="s1">&#39;plane strain&#39;</span><span class="p">:</span>
        <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="n">nu</span><span class="o">*</span><span class="p">(</span><span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">])</span>              <span class="c1"># zz</span>

    <span class="k">return</span> <span class="n">sigma</span></div>


<span class="k">class</span> <span class="nc">IrwinStressField</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Calculator to return Irwin near-tip stress field at atomic sites</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">K</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">x0</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">y0</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">sxx0</span><span class="o">=</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">syy0</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">sxy0</span><span class="o">=</span><span class="mf">0.</span><span class="p">,</span> <span class="n">nu</span><span class="o">=</span><span class="mf">0.5</span><span class="p">,</span>
                 <span class="n">stress_state</span><span class="o">=</span><span class="s1">&#39;plane strain&#39;</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">K</span> <span class="o">=</span> <span class="n">K</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">x0</span> <span class="o">=</span> <span class="n">x0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">y0</span> <span class="o">=</span> <span class="n">y0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">sxx0</span> <span class="o">=</span> <span class="n">sxx0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">syy0</span> <span class="o">=</span> <span class="n">syy0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">sxy0</span> <span class="o">=</span> <span class="n">sxy0</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">nu</span> <span class="o">=</span> <span class="n">nu</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">stress_state</span> <span class="o">=</span> <span class="n">stress_state</span>

    <span class="k">def</span> <span class="nf">get_stresses</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">atoms</span><span class="p">):</span>
        <span class="n">K</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">K</span>
        <span class="k">if</span> <span class="n">K</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">K</span> <span class="o">=</span> <span class="n">get_stress_intensity_factor</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>

        <span class="n">x0</span><span class="p">,</span> <span class="n">y0</span> <span class="o">=</span> <span class="bp">self</span><span class="o">.</span><span class="n">x0</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">y0</span>
        <span class="k">if</span> <span class="n">x0</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">x0</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">if</span> <span class="n">y0</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">y0</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span>

        <span class="n">x</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span>
        <span class="n">y</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span>
        <span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">x</span> <span class="o">-</span> <span class="n">x0</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="p">(</span><span class="n">y</span> <span class="o">-</span> <span class="n">y0</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
        <span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arctan2</span><span class="p">(</span><span class="n">y</span> <span class="o">-</span> <span class="n">y0</span><span class="p">,</span> <span class="n">x</span> <span class="o">-</span> <span class="n">x0</span><span class="p">)</span>

        <span class="n">sigma</span> <span class="o">=</span> <span class="n">irwin_modeI_crack_tip_stress_field</span><span class="p">(</span><span class="n">K</span><span class="p">,</span> <span class="n">r</span><span class="p">,</span> <span class="n">t</span><span class="p">,</span> <span class="bp">self</span><span class="o">.</span><span class="n">nu</span><span class="p">,</span>
                                                   <span class="bp">self</span><span class="o">.</span><span class="n">stress_state</span><span class="p">)</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+=</span> <span class="bp">self</span><span class="o">.</span><span class="n">sxx0</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="bp">self</span><span class="o">.</span><span class="n">syy0</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="bp">self</span><span class="o">.</span><span class="n">sxy0</span>
        <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+=</span> <span class="bp">self</span><span class="o">.</span><span class="n">sxy0</span>

        <span class="k">return</span> <span class="n">sigma</span>


<div class="viewcode-block" id="strain_to_G"><a class="viewcode-back" href="../../crack.html#quippy.crack.strain_to_G">[docs]</a><span class="k">def</span> <span class="nf">strain_to_G</span><span class="p">(</span><span class="n">strain</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="n">nu</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert from strain to energy release rate G for thin strip geometry</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    strain : float</span>
<span class="sd">       Dimensionless ratio ``(current_height - orig_height)/orig_height``</span>
<span class="sd">    E : float</span>
<span class="sd">       Young&#39;s modulus relevant for a pull in y direction sigma_yy/eps_yy</span>
<span class="sd">    nu : float</span>
<span class="sd">       Poission ratio -eps_yy/eps_xx</span>
<span class="sd">    orig_height : float</span>
<span class="sd">       Unstrained height of slab</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    G : float</span>
<span class="sd">       Energy release rate in units consistent with input</span>
<span class="sd">       (i.e. in eV/A**2 if eV/A/fs units used)</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="n">E</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.0</span> <span class="o">-</span> <span class="n">nu</span> <span class="o">*</span> <span class="n">nu</span><span class="p">)</span> <span class="o">*</span> <span class="n">strain</span> <span class="o">*</span> <span class="n">strain</span> <span class="o">*</span> <span class="n">orig_height</span></div>


<div class="viewcode-block" id="G_to_strain"><a class="viewcode-back" href="../../crack.html#quippy.crack.G_to_strain">[docs]</a><span class="k">def</span> <span class="nf">G_to_strain</span><span class="p">(</span><span class="n">G</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="n">nu</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Convert from energy release rate G to strain for thin strip geometry</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    G : float</span>
<span class="sd">       Energy release rate in units consistent with `E` and `orig_height`</span>
<span class="sd">    E : float</span>
<span class="sd">       Young&#39;s modulus relevant for a pull in y direction sigma_yy/eps_yy</span>
<span class="sd">    nu : float</span>
<span class="sd">       Poission ratio -eps_yy/eps_xx</span>
<span class="sd">    orig_height : float</span>
<span class="sd">       Unstrained height of slab</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    strain : float</span>
<span class="sd">       Dimensionless ratio ``(current_height - orig_height)/orig_height``</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="k">return</span> <span class="n">sqrt</span><span class="p">(</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">G</span> <span class="o">*</span> <span class="p">(</span><span class="mf">1.0</span> <span class="o">-</span> <span class="n">nu</span> <span class="o">*</span> <span class="n">nu</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">E</span> <span class="o">*</span> <span class="n">orig_height</span><span class="p">))</span></div>


<div class="viewcode-block" id="get_strain"><a class="viewcode-back" href="../../crack.html#quippy.crack.get_strain">[docs]</a><span class="k">def</span> <span class="nf">get_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return the current strain on thin strip configuration `atoms`</span>

<span class="sd">    Requires unstrained height of slab to be stored as ``OrigHeight``</span>
<span class="sd">    key in ``atoms.info`` dictionary.</span>

<span class="sd">    Also updates value stored in ``atoms.info``.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">orig_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>
    <span class="n">current_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
    <span class="n">strain</span> <span class="o">=</span> <span class="n">current_height</span> <span class="o">/</span> <span class="n">orig_height</span> <span class="o">-</span> <span class="mf">1.0</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">strain</span>
    <span class="k">return</span> <span class="n">strain</span></div>


<div class="viewcode-block" id="get_energy_release_rate"><a class="viewcode-back" href="../../crack.html#quippy.crack.get_energy_release_rate">[docs]</a><span class="k">def</span> <span class="nf">get_energy_release_rate</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return the current energy release rate G for `atoms</span>

<span class="sd">    Also updates value stored in ``atoms.info`` dictionary.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">current_strain</span> <span class="o">=</span> <span class="n">get_strain</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
    <span class="n">orig_height</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;OrigHeight&#39;</span><span class="p">]</span>
    <span class="n">E</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;YoungsModulus&#39;</span><span class="p">]</span>
    <span class="n">nu</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yx&#39;</span><span class="p">]</span>
    <span class="n">G</span> <span class="o">=</span> <span class="n">strain_to_G</span><span class="p">(</span><span class="n">current_strain</span><span class="p">,</span> <span class="n">E</span><span class="p">,</span> <span class="n">nu</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;G&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">G</span>
    <span class="k">return</span> <span class="n">G</span></div>


<div class="viewcode-block" id="get_stress_intensity_factor"><a class="viewcode-back" href="../../crack.html#quippy.crack.get_stress_intensity_factor">[docs]</a><span class="k">def</span> <span class="nf">get_stress_intensity_factor</span><span class="p">(</span><span class="n">atoms</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return stress intensity factor K_I</span>

<span class="sd">    Also updates value stored in ``atoms.info`` dictionary.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">G</span> <span class="o">=</span> <span class="n">get_energy_release_rate</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>

    <span class="n">E</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;YoungsModulus&#39;</span><span class="p">]</span>
    <span class="n">nu</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yx&#39;</span><span class="p">]</span>

    <span class="n">Ep</span> <span class="o">=</span> <span class="n">E</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">nu</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
    <span class="n">K</span> <span class="o">=</span> <span class="n">sqrt</span><span class="p">(</span><span class="n">G</span><span class="o">/</span><span class="n">Ep</span><span class="p">)</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">K</span>
    <span class="k">return</span> <span class="n">K</span></div>


<div class="viewcode-block" id="fit_crack_stress_field"><a class="viewcode-back" href="../../crack.html#quippy.crack.fit_crack_stress_field">[docs]</a><span class="k">def</span> <span class="nf">fit_crack_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">r_range</span><span class="o">=</span><span class="p">(</span><span class="mf">0.</span><span class="p">,</span> <span class="mf">50.</span><span class="p">),</span> <span class="n">initial_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">fix_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                           <span class="n">sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_decay</span><span class="o">=</span><span class="mf">0.005</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">verbose</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Perform a least squares fit of near-tip stress field to Irwin solution</span>

<span class="sd">    Stresses on the atoms are fit to the Irwin K-field singular crack tip</span>
<span class="sd">    solution, allowingthe crack position, stress intensity factor and</span>
<span class="sd">    far-field stress components to vary during the fit.</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    atoms : :class:`~.Atoms` object</span>
<span class="sd">       Crack system. For the initial fit, the following keys are used</span>
<span class="sd">       from the :attr:`~Atoms.info` dictionary:</span>

<span class="sd">          - ``YoungsModulus``</span>
<span class="sd">          - ``PossionRatio_yx``</span>
<span class="sd">          - ``G`` --- current energy release rate</span>
<span class="sd">          - ``strain`` --- current applied strain</span>
<span class="sd">          - ``CrackPos`` --- initial guess for crack tip position</span>

<span class="sd">       The initial guesses for the stress intensity factor ``K`` are</span>
<span class="sd">       far-field stress ``sigma0`` are computed from</span>
<span class="sd">       ``YoungsModulus``, ``PoissonRatio_yx``, ``G`` and ``strain``,</span>
<span class="sd">       assuming plane strain in thin strip boundary conditions.</span>

<span class="sd">       On exit, new ``K``, ``sigma0`` and ``CrackPos`` entries are set</span>
<span class="sd">       in the :attr:`~Atoms.info` dictionary. These values are then</span>
<span class="sd">       used as starting guesses for subsequent fits.</span>

<span class="sd">    r_range : sequence of two floats, optional</span>
<span class="sd">       If present, restrict the stress fit to an annular region</span>
<span class="sd">       ``r_range[0] &lt;= r &lt; r_range[1]``, centred on the previous crack</span>
<span class="sd">       position (from the ``CrackPos`` entry in ``atoms.info``). If</span>
<span class="sd">       r_range is ``None``, fit is carried out for all atoms.</span>

<span class="sd">    initial_params : dict</span>
<span class="sd">       Names and initial values of parameters. Missing initial values</span>
<span class="sd">       are guessed from Atoms object.</span>

<span class="sd">    fix_params : dict</span>
<span class="sd">       Names and values of parameters to fix during the fit,</span>
<span class="sd">       e.g. ``{y0: 0.0}`` to constrain the fit to the line y=0</span>

<span class="sd">    sigma : None or array with shape (len(atoms), 3, 3)</span>
<span class="sd">       Explicitly provide the per-atom stresses. Avoids calling Atoms&#39;</span>
<span class="sd">       calculators :meth:`~.get_stresses` method.</span>

<span class="sd">    avg_sigma : None or array with shape (len(atoms), 3, 3)</span>
<span class="sd">       If present, use this array to accumulate the time-averaged</span>
<span class="sd">       stress field. Useful when processing a trajectory.</span>

<span class="sd">    avg_decay : real</span>
<span class="sd">       Factor by which average stress is attenuated at each step.</span>
<span class="sd">       Should be set to ``dt/tau`` where ``dt`` is MD time-step</span>
<span class="sd">       and ``tau`` is a characteristic averaging time.</span>

<span class="sd">    calc : Calculator object, optional</span>
<span class="sd">       If present, override the calculator used to compute stresses</span>
<span class="sd">       on the atoms. Default is ``atoms.get_calculator``.</span>

<span class="sd">       To use the atom resolved stress tensor pass an instance of the</span>
<span class="sd">       :class:`~quippy.elasticity.AtomResolvedStressField` class.</span>

<span class="sd">    verbose : bool, optional</span>
<span class="sd">       If set to True, print additional information about the fit.</span>

<span class="sd">    Returns</span>
<span class="sd">    -------</span>
<span class="sd">    params : dict with keys ``[K, x0, y0, sxx0, syy0, sxy0]``</span>
<span class="sd">       Fitted parameters, in a form suitable for passin</span>
<span class="sd">       :class:`IrwinStressField` constructor. These are the stress intensity</span>
<span class="sd">       factor `K`, the centre of the stress field ``(x0, y0)``, and the</span>
<span class="sd">       far field contribution to the stress ``(sxx0, syy0, sxy0)``.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">params</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="k">if</span> <span class="n">initial_params</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">params</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">initial_params</span><span class="p">)</span>

    <span class="k">if</span> <span class="s1">&#39;K&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span><span class="p">:</span>
       <span class="c1"># Guess for stress intensity factor K</span>
       <span class="k">if</span> <span class="s1">&#39;K&#39;</span> <span class="ow">in</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">:</span>
           <span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span>
       <span class="k">else</span><span class="p">:</span>
           <span class="k">try</span><span class="p">:</span>
               <span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">get_stress_intensity_factor</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
           <span class="k">except</span> <span class="ne">KeyError</span><span class="p">:</span>
               <span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">1.0</span><span class="o">*</span><span class="n">MPA_SQRT_M</span>

    <span class="k">if</span> <span class="s1">&#39;sxx0&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span> <span class="ow">or</span> <span class="s1">&#39;syy0&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span> <span class="ow">or</span> <span class="s1">&#39;sxy0&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span><span class="p">:</span>
       <span class="c1"># Guess for far-field stress</span>
       <span class="k">if</span> <span class="s1">&#39;sigma0&#39;</span> <span class="ow">in</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">:</span>
          <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;sigma0&#39;</span><span class="p">]</span>
       <span class="k">else</span><span class="p">:</span>
          <span class="k">try</span><span class="p">:</span>
              <span class="n">E</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;YoungsModulus&#39;</span><span class="p">]</span>
              <span class="n">nu</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;PoissonRatio_yx&#39;</span><span class="p">]</span>
              <span class="n">Ep</span> <span class="o">=</span> <span class="n">E</span><span class="o">/</span><span class="p">(</span><span class="mi">1</span><span class="o">-</span><span class="n">nu</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">Ep</span><span class="o">*</span><span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;strain&#39;</span><span class="p">]</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">nu</span><span class="o">*</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">]</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
          <span class="k">except</span> <span class="ne">KeyError</span><span class="p">:</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>

    <span class="k">if</span> <span class="s1">&#39;x0&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span> <span class="ow">or</span> <span class="s1">&#39;y0&#39;</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">params</span><span class="p">:</span>
       <span class="c1"># Guess for crack position</span>
       <span class="k">try</span><span class="p">:</span>
           <span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">],</span> <span class="n">_</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span>
       <span class="k">except</span> <span class="ne">KeyError</span><span class="p">:</span>
           <span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span> <span class="o">+</span>
                           <span class="p">(</span><span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">())</span><span class="o">/</span><span class="mf">3.0</span><span class="p">)</span>
           <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>

    <span class="c1"># Override any fixed parameters</span>
    <span class="k">if</span> <span class="n">fix_params</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">fix_params</span> <span class="o">=</span> <span class="p">{}</span>
    <span class="n">params</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">fix_params</span><span class="p">)</span>

    <span class="n">x</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">x</span> <span class="o">-</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="p">(</span><span class="n">y</span> <span class="o">-</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">])</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>

    <span class="c1"># Get local stresses</span>
    <span class="k">if</span> <span class="n">sigma</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="k">if</span> <span class="n">calc</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
           <span class="n">calc</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">get_calculator</span><span class="p">()</span>
       <span class="n">sigma</span> <span class="o">=</span> <span class="n">calc</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>

    <span class="c1"># Update avg_sigma in place</span>
    <span class="k">if</span> <span class="n">avg_sigma</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">avg_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">avg_decay</span><span class="p">)</span><span class="o">*</span><span class="n">avg_sigma</span> <span class="o">+</span> <span class="p">(</span><span class="mf">1.0</span> <span class="o">-</span> <span class="n">np</span><span class="o">.</span><span class="n">exp</span><span class="p">(</span><span class="o">-</span><span class="n">avg_decay</span><span class="p">))</span><span class="o">*</span><span class="n">sigma</span>
       <span class="n">sigma</span> <span class="o">=</span> <span class="n">avg_sigma</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>

    <span class="c1"># Zero components out of the xy plane</span>
    <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
    <span class="n">sigma</span><span class="p">[:,</span><span class="mi">0</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
    <span class="n">sigma</span><span class="p">[:,</span><span class="mi">2</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
    <span class="n">sigma</span><span class="p">[:,</span><span class="mi">1</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>
    <span class="n">sigma</span><span class="p">[:,</span><span class="mf">2.1</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.0</span>

    <span class="n">mask</span> <span class="o">=</span> <span class="bp">Ellipsis</span> <span class="c1"># all atoms</span>
    <span class="k">if</span> <span class="n">r_range</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
        <span class="n">rmin</span><span class="p">,</span> <span class="n">rmax</span> <span class="o">=</span> <span class="n">r_range</span>
        <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">r</span> <span class="o">&gt;</span> <span class="n">rmin</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">r</span> <span class="o">&lt;</span> <span class="n">rmax</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
       <span class="nb">print</span> <span class="s1">&#39;Fitting on </span><span class="si">%r</span><span class="s1"> atoms&#39;</span> <span class="o">%</span> <span class="n">sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">shape</span>

    <span class="k">def</span> <span class="nf">objective_function</span><span class="p">(</span><span class="n">params</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">var_params</span><span class="p">):</span>
        <span class="n">params</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span><span class="nb">zip</span><span class="p">(</span><span class="n">var_params</span><span class="p">,</span> <span class="n">params</span><span class="p">))</span>
        <span class="k">if</span> <span class="n">fix_params</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">params</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">fix_params</span><span class="p">)</span>
        <span class="n">irwin_sigma</span> <span class="o">=</span> <span class="n">IrwinStressField</span><span class="p">(</span><span class="o">**</span><span class="n">params</span><span class="p">)</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">(</span><span class="n">atoms</span><span class="p">)</span>
        <span class="n">delta_sigma</span> <span class="o">=</span> <span class="n">sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,:,:]</span> <span class="o">-</span> <span class="n">irwin_sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,:,:]</span>
        <span class="k">return</span> <span class="n">delta_sigma</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="n">delta_sigma</span><span class="o">.</span><span class="n">size</span><span class="p">)</span>

    <span class="c1"># names and values of parameters which can vary in this fit</span>
    <span class="n">var_params</span> <span class="o">=</span> <span class="nb">sorted</span><span class="p">([</span><span class="n">key</span> <span class="k">for</span> <span class="n">key</span> <span class="ow">in</span> <span class="n">params</span><span class="o">.</span><span class="n">keys</span><span class="p">()</span> <span class="k">if</span> <span class="n">key</span> <span class="ow">not</span> <span class="ow">in</span> <span class="n">fix_params</span><span class="o">.</span><span class="n">keys</span><span class="p">()</span> <span class="p">])</span>
    <span class="n">initial_params</span> <span class="o">=</span> <span class="p">[</span><span class="n">params</span><span class="p">[</span><span class="n">key</span><span class="p">]</span> <span class="k">for</span> <span class="n">key</span> <span class="ow">in</span> <span class="n">var_params</span><span class="p">]</span>

    <span class="kn">from</span> <span class="nn">scipy.optimize</span> <span class="k">import</span> <span class="n">leastsq</span>
    <span class="n">fitted_params</span><span class="p">,</span> <span class="n">cov</span><span class="p">,</span> <span class="n">infodict</span><span class="p">,</span> <span class="n">mesg</span><span class="p">,</span> <span class="n">success</span> <span class="o">=</span> <span class="n">leastsq</span><span class="p">(</span><span class="n">objective_function</span><span class="p">,</span>
                                                         <span class="n">initial_params</span><span class="p">,</span>
                                                         <span class="n">args</span><span class="o">=</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">var_params</span><span class="p">),</span>
                                                         <span class="n">full_output</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>

    <span class="n">params</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span><span class="nb">zip</span><span class="p">(</span><span class="n">var_params</span><span class="p">,</span> <span class="n">fitted_params</span><span class="p">))</span>
    <span class="n">params</span><span class="o">.</span><span class="n">update</span><span class="p">(</span><span class="n">fix_params</span><span class="p">)</span>

    <span class="c1"># estimate variance in parameter estimates</span>
    <span class="k">if</span> <span class="n">cov</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="c1"># singular covariance matrix</span>
       <span class="n">err</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span><span class="nb">zip</span><span class="p">(</span><span class="n">var_params</span><span class="p">,</span> <span class="p">[</span><span class="mf">0.</span><span class="p">]</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">fitted_params</span><span class="p">)))</span>
    <span class="k">else</span><span class="p">:</span>
       <span class="n">s_sq</span> <span class="o">=</span> <span class="p">(</span><span class="n">objective_function</span><span class="p">(</span><span class="n">fitted_params</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span> <span class="n">var_params</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span><span class="o">.</span><span class="n">sum</span><span class="p">()</span><span class="o">/</span><span class="p">(</span><span class="n">sigma</span><span class="o">.</span><span class="n">size</span><span class="o">-</span><span class="nb">len</span><span class="p">(</span><span class="n">fitted_params</span><span class="p">))</span>
       <span class="n">cov</span> <span class="o">=</span> <span class="n">cov</span> <span class="o">*</span> <span class="n">s_sq</span>
       <span class="n">err</span> <span class="o">=</span> <span class="nb">dict</span><span class="p">(</span><span class="nb">zip</span><span class="p">(</span><span class="n">var_params</span><span class="p">,</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">diag</span><span class="p">(</span><span class="n">cov</span><span class="p">))))</span>

    <span class="k">if</span> <span class="n">verbose</span><span class="p">:</span>
       <span class="nb">print</span> <span class="s1">&#39;K = </span><span class="si">%.3f</span><span class="s1"> MPa sqrt(m)&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span><span class="o">/</span><span class="n">MPA_SQRT_M</span><span class="p">)</span>
       <span class="nb">print</span> <span class="s1">&#39;sigma^0_{xx,yy,xy} = (</span><span class="si">%.1f</span><span class="s1">, </span><span class="si">%.1f</span><span class="s1">, </span><span class="si">%.1f</span><span class="s1">) GPa&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">]</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">,</span>
                                                              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">]</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">,</span>
                                                              <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">]</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">)</span>
       <span class="nb">print</span> <span class="s1">&#39;Crack position (x0, y0) = (</span><span class="si">%.1f</span><span class="s1">, </span><span class="si">%.1f</span><span class="s1">) A&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">])</span>

    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">]</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;sigma0&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">])</span>
    <span class="n">atoms</span><span class="o">.</span><span class="n">info</span><span class="p">[</span><span class="s1">&#39;CrackPos&#39;</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">((</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">],</span> <span class="n">atoms</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">/</span><span class="mf">2.0</span><span class="p">))</span>

    <span class="k">return</span> <span class="n">params</span><span class="p">,</span> <span class="n">err</span></div>



<div class="viewcode-block" id="find_crack_tip_stress_field"><a class="viewcode-back" href="../../crack.html#quippy.crack.find_crack_tip_stress_field">[docs]</a><span class="k">def</span> <span class="nf">find_crack_tip_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">r_range</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">initial_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">fix_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                                <span class="n">sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_decay</span><span class="o">=</span><span class="mf">0.005</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Find the position of the crack tip by fitting to the Irwin `K`-field solution</span>

<span class="sd">    Fit is carried out using :func:`fit_crack_stress_field`, and parameters</span>
<span class="sd">    have the same meaning as there.</span>

<span class="sd">    See also</span>
<span class="sd">    --------</span>
<span class="sd">    fit_crack_stress_field</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">params</span><span class="p">,</span> <span class="n">err</span> <span class="o">=</span> <span class="n">fit_crack_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">r_range</span><span class="p">,</span> <span class="n">initial_params</span><span class="p">,</span> <span class="n">fix_params</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span>
                                         <span class="n">avg_sigma</span><span class="p">,</span> <span class="n">avg_decay</span><span class="p">,</span> <span class="n">calc</span><span class="p">)</span>

    <span class="k">return</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">((</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">],</span> <span class="n">atoms</span><span class="o">.</span><span class="n">cell</span><span class="p">[</span><span class="mi">2</span><span class="p">,</span><span class="mi">2</span><span class="p">]</span><span class="o">/</span><span class="mf">2.0</span><span class="p">))</span></div>


<div class="viewcode-block" id="plot_stress_fields"><a class="viewcode-back" href="../../crack.html#quippy.crack.plot_stress_fields">[docs]</a><span class="k">def</span> <span class="nf">plot_stress_fields</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">r_range</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">initial_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">fix_params</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span>
                       <span class="n">sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_sigma</span><span class="o">=</span><span class="kc">None</span><span class="p">,</span> <span class="n">avg_decay</span><span class="o">=</span><span class="mf">0.005</span><span class="p">,</span> <span class="n">calc</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Fit and plot atomistic and continuum stress fields</span>

<span class="sd">    Firstly a fit to the Irwin `K`-field solution is carried out using</span>
<span class="sd">    :func:`fit_crack_stress_field`, and parameters have the same</span>
<span class="sd">    meaning as for that function. Then plots of the</span>
<span class="sd">    :math:`\sigma_{xx}`, :math:`\sigma_{yy}`, :math:`\sigma_{xy}`</span>
<span class="sd">    fields are produced for atomistic and continuum cases, and for the</span>
<span class="sd">    residual error after fitting.</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="kn">from</span> <span class="nn">pylab</span> <span class="k">import</span> <span class="n">griddata</span><span class="p">,</span> <span class="n">meshgrid</span><span class="p">,</span> <span class="n">subplot</span><span class="p">,</span> <span class="n">cla</span><span class="p">,</span> <span class="n">contourf</span><span class="p">,</span> <span class="n">colorbar</span><span class="p">,</span> <span class="n">draw</span><span class="p">,</span> <span class="n">title</span><span class="p">,</span> <span class="n">clf</span><span class="p">,</span> <span class="n">gca</span>

    <span class="n">params</span><span class="p">,</span> <span class="n">err</span> <span class="o">=</span> <span class="n">fit_crack_stress_field</span><span class="p">(</span><span class="n">atoms</span><span class="p">,</span> <span class="n">r_range</span><span class="p">,</span> <span class="n">initial_params</span><span class="p">,</span> <span class="n">fix_params</span><span class="p">,</span> <span class="n">sigma</span><span class="p">,</span>
                                         <span class="n">avg_sigma</span><span class="p">,</span> <span class="n">avg_decay</span><span class="p">,</span> <span class="n">calc</span><span class="p">)</span>

    <span class="n">K</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">y0</span><span class="p">,</span> <span class="n">sxx0</span><span class="p">,</span> <span class="n">syy0</span><span class="p">,</span> <span class="n">sxy0</span> <span class="o">=</span> <span class="p">(</span><span class="n">params</span><span class="p">[</span><span class="s1">&#39;K&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;x0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;y0&#39;</span><span class="p">],</span>
                                   <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxx0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;syy0&#39;</span><span class="p">],</span> <span class="n">params</span><span class="p">[</span><span class="s1">&#39;sxy0&#39;</span><span class="p">])</span>

    <span class="n">x</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">0</span><span class="p">]</span>
    <span class="n">y</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">positions</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span>

    <span class="n">X</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">((</span><span class="n">x</span><span class="o">-</span><span class="n">x0</span><span class="p">)</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="p">(</span><span class="n">x</span><span class="o">-</span><span class="n">x0</span><span class="p">)</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span> <span class="mi">500</span><span class="p">)</span>
    <span class="n">Y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">((</span><span class="n">y</span><span class="o">-</span><span class="n">y0</span><span class="p">)</span><span class="o">.</span><span class="n">min</span><span class="p">(),</span> <span class="p">(</span><span class="n">y</span><span class="o">-</span><span class="n">y0</span><span class="p">)</span><span class="o">.</span><span class="n">max</span><span class="p">(),</span> <span class="mi">500</span><span class="p">)</span>

    <span class="n">t</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arctan2</span><span class="p">(</span><span class="n">y</span><span class="o">-</span><span class="n">y0</span><span class="p">,</span> <span class="n">x</span><span class="o">-</span><span class="n">x0</span><span class="p">)</span>
    <span class="n">r</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">((</span><span class="n">x</span><span class="o">-</span><span class="n">x0</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span> <span class="o">+</span> <span class="p">(</span><span class="n">y</span><span class="o">-</span><span class="n">y0</span><span class="p">)</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>

    <span class="k">if</span> <span class="n">r_range</span> <span class="ow">is</span> <span class="ow">not</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">rmin</span><span class="p">,</span> <span class="n">rmax</span> <span class="o">=</span> <span class="n">r_range</span>
       <span class="n">mask</span> <span class="o">=</span> <span class="p">(</span><span class="n">r</span> <span class="o">&gt;</span> <span class="n">rmin</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">r</span> <span class="o">&lt;</span> <span class="n">rmax</span><span class="p">)</span>
    <span class="k">else</span><span class="p">:</span>
       <span class="n">mask</span> <span class="o">=</span> <span class="bp">Ellipsis</span>

    <span class="n">atom_sigma</span> <span class="o">=</span> <span class="n">sigma</span>
    <span class="k">if</span> <span class="n">atom_sigma</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
       <span class="n">atom_sigma</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">get_stresses</span><span class="p">()</span>

    <span class="n">grid_sigma</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">dstack</span><span class="p">([</span><span class="n">griddata</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">x0</span><span class="p">,</span> <span class="n">y</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">y0</span><span class="p">,</span> <span class="n">atom_sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">],</span> <span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">),</span>
                            <span class="n">griddata</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">x0</span><span class="p">,</span> <span class="n">y</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">y0</span><span class="p">,</span> <span class="n">atom_sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">),</span>
                            <span class="n">griddata</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">x0</span><span class="p">,</span> <span class="n">y</span><span class="p">[</span><span class="n">mask</span><span class="p">]</span><span class="o">-</span><span class="n">y0</span><span class="p">,</span> <span class="n">atom_sigma</span><span class="p">[</span><span class="n">mask</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">],</span> <span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">)])</span>

    <span class="n">X</span><span class="p">,</span> <span class="n">Y</span> <span class="o">=</span> <span class="n">meshgrid</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">)</span>
    <span class="n">R</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">X</span><span class="o">**</span><span class="mi">2</span><span class="o">+</span><span class="n">Y</span><span class="o">**</span><span class="mi">2</span><span class="p">)</span>
    <span class="n">T</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">arctan2</span><span class="p">(</span><span class="n">Y</span><span class="p">,</span> <span class="n">X</span><span class="p">)</span>

    <span class="n">grid_sigma</span><span class="p">[((</span><span class="n">R</span> <span class="o">&lt;</span> <span class="n">rmin</span><span class="p">)</span> <span class="o">|</span> <span class="p">(</span><span class="n">R</span> <span class="o">&gt;</span> <span class="n">rmax</span><span class="p">)),:]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">nan</span> <span class="c1"># mask outside fitting region</span>

    <span class="n">irwin_sigma</span> <span class="o">=</span> <span class="n">irwin_modeI_crack_tip_stress_field</span><span class="p">(</span><span class="n">K</span><span class="p">,</span> <span class="n">R</span><span class="p">,</span> <span class="n">T</span><span class="p">,</span> <span class="n">x0</span><span class="p">,</span> <span class="n">y0</span><span class="p">)</span>
    <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+=</span> <span class="n">sxx0</span>
    <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="n">syy0</span>
    <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="n">sxy0</span>
    <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="mi">1</span><span class="p">,</span><span class="mi">0</span><span class="p">]</span> <span class="o">+=</span> <span class="n">sxy0</span>
    <span class="n">irwin_sigma</span> <span class="o">=</span> <span class="n">ma</span><span class="o">.</span><span class="n">masked_array</span><span class="p">(</span><span class="n">irwin_sigma</span><span class="p">,</span> <span class="n">mask</span><span class="o">=</span><span class="n">grid_sigma</span><span class="o">.</span><span class="n">mask</span><span class="p">)</span>

    <span class="n">irwin_sigma</span><span class="p">[((</span><span class="n">R</span> <span class="o">&lt;</span> <span class="n">rmin</span><span class="p">)</span> <span class="o">|</span> <span class="p">(</span><span class="n">R</span> <span class="o">&gt;</span> <span class="n">rmax</span><span class="p">)),:,:]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">nan</span> <span class="c1"># mask outside fitting region</span>

    <span class="n">contours</span> <span class="o">=</span> <span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">20</span><span class="p">,</span> <span class="mi">10</span><span class="p">),</span>
                <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">20</span><span class="p">,</span> <span class="mi">10</span><span class="p">),</span>
                <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">10</span><span class="p">,</span><span class="mi">10</span><span class="p">,</span> <span class="mi">10</span><span class="p">)]</span>

    <span class="n">dcontours</span> <span class="o">=</span> <span class="p">[</span><span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">10</span><span class="p">),</span>
                <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">10</span><span class="p">),</span>
                <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">5</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">10</span><span class="p">)]</span>

    <span class="n">clf</span><span class="p">()</span>
    <span class="k">for</span> <span class="n">i</span><span class="p">,</span> <span class="p">(</span><span class="n">ii</span><span class="p">,</span> <span class="n">jj</span><span class="p">),</span> <span class="n">label</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">3</span><span class="p">),</span>
                                  <span class="p">[(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">),</span> <span class="p">(</span><span class="mi">1</span><span class="p">,</span><span class="mi">1</span><span class="p">),</span> <span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">1</span><span class="p">)],</span>
                                  <span class="p">[</span><span class="s1">&#39;\sigma_</span><span class="si">{xx}</span><span class="s1">&#39;</span><span class="p">,</span> <span class="sa">r</span><span class="s1">&#39;\sigma_</span><span class="si">{yy}</span><span class="s1">&#39;</span><span class="p">,</span> <span class="sa">r</span><span class="s1">&#39;\sigma_</span><span class="si">{xy}</span><span class="s1">&#39;</span><span class="p">]):</span>
        <span class="n">subplot</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">)</span>
        <span class="n">gca</span><span class="p">()</span><span class="o">.</span><span class="n">set_aspect</span><span class="p">(</span><span class="s1">&#39;equal&#39;</span><span class="p">)</span>
        <span class="n">contourf</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">,</span> <span class="n">grid_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">i</span><span class="p">]</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">,</span> <span class="n">contours</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
        <span class="n">colorbar</span><span class="p">()</span>
        <span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;$</span><span class="si">%s</span><span class="s1">^\mathrm</span><span class="si">{atom}</span><span class="s1">$&#39;</span> <span class="o">%</span> <span class="n">label</span><span class="p">)</span>
        <span class="n">draw</span><span class="p">()</span>

        <span class="n">subplot</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="n">i</span><span class="o">+</span><span class="mi">4</span><span class="p">)</span>
        <span class="n">gca</span><span class="p">()</span><span class="o">.</span><span class="n">set_aspect</span><span class="p">(</span><span class="s1">&#39;equal&#39;</span><span class="p">)</span>
        <span class="n">contourf</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">,</span> <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">ii</span><span class="p">,</span><span class="n">jj</span><span class="p">]</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">,</span> <span class="n">contours</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
        <span class="n">colorbar</span><span class="p">()</span>
        <span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;$</span><span class="si">%s</span><span class="s1">^\mathrm</span><span class="si">{Irwin}</span><span class="s1">$&#39;</span> <span class="o">%</span> <span class="n">label</span><span class="p">)</span>
        <span class="n">draw</span><span class="p">()</span>

        <span class="n">subplot</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">3</span><span class="p">,</span><span class="n">i</span><span class="o">+</span><span class="mi">7</span><span class="p">)</span>
        <span class="n">gca</span><span class="p">()</span><span class="o">.</span><span class="n">set_aspect</span><span class="p">(</span><span class="s1">&#39;equal&#39;</span><span class="p">)</span>
        <span class="n">contourf</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">,</span> <span class="nb">abs</span><span class="p">(</span><span class="n">grid_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span>
                           <span class="n">irwin_sigma</span><span class="p">[</span><span class="o">...</span><span class="p">,</span><span class="n">ii</span><span class="p">,</span><span class="n">jj</span><span class="p">])</span><span class="o">*</span><span class="n">EV_A3_IN_GPA</span><span class="p">,</span> <span class="n">dcontours</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
        <span class="n">colorbar</span><span class="p">()</span>
        <span class="n">title</span><span class="p">(</span><span class="sa">r</span><span class="s1">&#39;$|</span><span class="si">%s</span><span class="s1">^\mathrm</span><span class="si">{atom}</span><span class="s1"> - </span><span class="si">%s</span><span class="s1">^\mathrm</span><span class="si">{Irwin}</span><span class="s1">|$&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">label</span><span class="p">,</span> <span class="n">label</span><span class="p">))</span>
        <span class="n">draw</span><span class="p">()</span></div>



<div class="viewcode-block" id="thin_strip_displacement_y"><a class="viewcode-back" href="../../crack.html#quippy.crack.thin_strip_displacement_y">[docs]</a><span class="k">def</span> <span class="nf">thin_strip_displacement_y</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">strain</span><span class="p">,</span> <span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Return vertical displacement ramp used to apply initial strain to slab</span>

<span class="sd">    Strain is increased from 0 to strain over distance :math:`a &lt;= x &lt;= b`.</span>
<span class="sd">    Region :math:`x &lt; a` is rigidly shifted up/down by ``strain*height/2``.</span>

<span class="sd">    Here is an example of how to use this function on an artificial</span>
<span class="sd">    2D square atomic lattice. The positions are plotted before (left)</span>
<span class="sd">    and after (right) applying the displacement, and the horizontal and</span>
<span class="sd">    vertical lines show the `strain` (red), `a` (green) and `b` (blue)</span>
<span class="sd">    parameters. ::</span>

<span class="sd">       import matplotlib.pyplot as plt</span>
<span class="sd">       import numpy as np</span>

<span class="sd">       w = 1; h = 1; strain = 0.1; a = -0.5; b =  0.0</span>
<span class="sd">       x = np.linspace(-w, w, 20)</span>
<span class="sd">       y = np.linspace(-h, h, 20)</span>
<span class="sd">       X, Y = np.meshgrid(x, y)</span>
<span class="sd">       u_y = thin_strip_displacement_y(X, Y, strain, a, b)</span>

<span class="sd">       for i, disp in enumerate([0, u_y]):</span>
<span class="sd">           plt.subplot(1,2,i+1)</span>
<span class="sd">           plt.scatter(X, Y + disp, c=&#39;k&#39;, s=5)</span>
<span class="sd">           for y in [-h, h]:</span>
<span class="sd">               plt.axhline(y, color=&#39;r&#39;, linewidth=2, linestyle=&#39;dashed&#39;)</span>
<span class="sd">               plt.axhline(y*(1+strain), color=&#39;r&#39;, linewidth=2)</span>
<span class="sd">           for x, c in zip([a, b], [&#39;g&#39;, &#39;b&#39;]):</span>
<span class="sd">               plt.axvline(x, color=c, linewidth=2)</span>

<span class="sd">    .. image:: thin-strip-displacement-y.png</span>
<span class="sd">       :width: 600</span>
<span class="sd">       :align: center</span>

<span class="sd">    Parameters</span>
<span class="sd">    ----------</span>
<span class="sd">    x : array</span>
<span class="sd">    y : array</span>
<span class="sd">       Atomic positions in unstrained slab, centered on origin x=0,y=0</span>
<span class="sd">    strain : float</span>
<span class="sd">       Far field strain to apply</span>
<span class="sd">    a : float</span>
<span class="sd">       x coordinate for beginning of strain ramp</span>
<span class="sd">    b : float</span>
<span class="sd">       x coordinate for end of strain ramp</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="n">u_y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros_like</span><span class="p">(</span><span class="n">y</span><span class="p">)</span>
    <span class="n">height</span> <span class="o">=</span> <span class="n">y</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">y</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>     <span class="c1"># measure height of slab</span>
    <span class="n">shift</span> <span class="o">=</span> <span class="n">strain</span> <span class="o">*</span> <span class="n">height</span> <span class="o">/</span> <span class="mf">2.0</span>  <span class="c1"># far behind crack, shift = strain*height/2</span>

    <span class="n">u_y</span><span class="p">[</span><span class="n">x</span> <span class="o">&lt;</span> <span class="n">a</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="n">x</span> <span class="o">&lt;</span> <span class="n">a</span><span class="p">])</span> <span class="o">*</span> <span class="n">shift</span>  <span class="c1"># region shift for x &lt; a</span>
    <span class="n">u_y</span><span class="p">[</span><span class="n">x</span> <span class="o">&gt;</span> <span class="n">b</span><span class="p">]</span> <span class="o">=</span> <span class="n">strain</span> <span class="o">*</span> <span class="n">y</span><span class="p">[</span><span class="n">x</span> <span class="o">&gt;</span> <span class="n">b</span><span class="p">]</span>          <span class="c1"># constant strain for x &gt; b</span>

    <span class="n">middle</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span> <span class="o">&gt;=</span> <span class="n">a</span><span class="p">)</span> <span class="o">&amp;</span> <span class="p">(</span><span class="n">x</span> <span class="o">&lt;=</span> <span class="n">b</span><span class="p">)</span>            <span class="c1"># interpolate for a &lt;= x &lt;= b</span>
    <span class="n">f</span> <span class="o">=</span> <span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">middle</span><span class="p">]</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">b</span> <span class="o">-</span> <span class="n">a</span><span class="p">)</span>
    <span class="n">u_y</span><span class="p">[</span><span class="n">middle</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span><span class="n">f</span> <span class="o">*</span> <span class="n">strain</span> <span class="o">*</span> <span class="n">y</span><span class="p">[</span><span class="n">middle</span><span class="p">]</span> <span class="o">+</span>
                   <span class="p">(</span><span class="mi">1</span> <span class="o">-</span> <span class="n">f</span><span class="p">)</span> <span class="o">*</span> <span class="n">shift</span> <span class="o">*</span> <span class="n">np</span><span class="o">.</span><span class="n">sign</span><span class="p">(</span><span class="n">y</span><span class="p">[</span><span class="n">middle</span><span class="p">]))</span>

    <span class="k">return</span> <span class="n">u_y</span></div>


<div class="viewcode-block" id="print_crack_system"><a class="viewcode-back" href="../../crack.html#quippy.crack.print_crack_system">[docs]</a><span class="k">def</span> <span class="nf">print_crack_system</span><span class="p">(</span><span class="n">crack_direction</span><span class="p">,</span> <span class="n">cleavage_plane</span><span class="p">,</span> <span class="n">crack_front</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Pretty printing of crack crystallographic coordinate system</span>

<span class="sd">    Specified by Miller indices for crack_direction (x),</span>
<span class="sd">    cleavage_plane (y) and crack_front (z), each of which should be</span>
<span class="sd">    a sequence of three floats</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">crack_direction</span> <span class="o">=</span> <span class="n">MillerDirection</span><span class="p">(</span><span class="n">crack_direction</span><span class="p">)</span>
    <span class="n">cleavage_plane</span> <span class="o">=</span> <span class="n">MillerPlane</span><span class="p">(</span><span class="n">cleavage_plane</span><span class="p">)</span>
    <span class="n">crack_front</span> <span class="o">=</span> <span class="n">MillerDirection</span><span class="p">(</span><span class="n">crack_front</span><span class="p">)</span>

    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Crack system              </span><span class="si">%s%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">cleavage_plane</span><span class="p">,</span> <span class="n">crack_front</span><span class="p">))</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Crack direction (x-axis)  </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">crack_direction</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Cleavage plane  (y-axis)  </span><span class="si">%s</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">cleavage_plane</span><span class="p">)</span>
    <span class="nb">print</span><span class="p">(</span><span class="s1">&#39;Crack front     (z-axis)  </span><span class="si">%s</span><span class="se">\n</span><span class="s1">&#39;</span> <span class="o">%</span> <span class="n">crack_front</span><span class="p">)</span></div>



<div class="viewcode-block" id="ConstantStrainRate"><a class="viewcode-back" href="../../crack.html#quippy.crack.ConstantStrainRate">[docs]</a><span class="k">class</span> <span class="nc">ConstantStrainRate</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Constraint which increments epsilon_yy at a constant strain rate</span>

<span class="sd">    Rescaling is applied only to atoms where `mask` is True (default is all atoms)</span>
<span class="sd">    &quot;&quot;&quot;</span>

    <span class="k">def</span> <span class="nf">__init__</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">orig_height</span><span class="p">,</span> <span class="n">delta_strain</span><span class="p">,</span> <span class="n">mask</span><span class="o">=</span><span class="kc">None</span><span class="p">):</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">orig_height</span> <span class="o">=</span> <span class="n">orig_height</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">delta_strain</span> <span class="o">=</span> <span class="n">delta_strain</span>
        <span class="k">if</span> <span class="n">mask</span> <span class="ow">is</span> <span class="kc">None</span><span class="p">:</span>
            <span class="n">mask</span> <span class="o">=</span> <span class="bp">Ellipsis</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">mask</span> <span class="o">=</span> <span class="n">mask</span>

    <span class="k">def</span> <span class="nf">adjust_forces</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">positions</span><span class="p">,</span> <span class="n">forces</span><span class="p">):</span>
        <span class="k">pass</span>

    <span class="k">def</span> <span class="nf">adjust_positions</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">newpos</span><span class="p">):</span>
        <span class="n">current_height</span> <span class="o">=</span> <span class="n">newpos</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">max</span><span class="p">()</span> <span class="o">-</span> <span class="n">newpos</span><span class="p">[:,</span> <span class="mi">1</span><span class="p">]</span><span class="o">.</span><span class="n">min</span><span class="p">()</span>
        <span class="n">current_strain</span> <span class="o">=</span> <span class="n">current_height</span> <span class="o">/</span> <span class="bp">self</span><span class="o">.</span><span class="n">orig_height</span> <span class="o">-</span> <span class="mf">1.0</span>
        <span class="n">new_strain</span> <span class="o">=</span> <span class="n">current_strain</span> <span class="o">+</span> <span class="bp">self</span><span class="o">.</span><span class="n">delta_strain</span>
        <span class="n">alpha</span> <span class="o">=</span> <span class="p">(</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">new_strain</span><span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="mf">1.0</span> <span class="o">+</span> <span class="n">current_strain</span><span class="p">)</span>
        <span class="n">newpos</span><span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">mask</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span> <span class="o">=</span> <span class="n">newpos</span><span class="p">[</span><span class="bp">self</span><span class="o">.</span><span class="n">mask</span><span class="p">,</span> <span class="mi">1</span><span class="p">]</span><span class="o">*</span><span class="n">alpha</span>

    <span class="k">def</span> <span class="nf">copy</span><span class="p">(</span><span class="bp">self</span><span class="p">):</span>
        <span class="k">return</span> <span class="n">ConstantStrainRate</span><span class="p">(</span><span class="bp">self</span><span class="o">.</span><span class="n">orig_height</span><span class="p">,</span>
                                  <span class="bp">self</span><span class="o">.</span><span class="n">delta_strain</span><span class="p">,</span>
                                  <span class="bp">self</span><span class="o">.</span><span class="n">mask</span><span class="p">)</span>

<div class="viewcode-block" id="ConstantStrainRate.apply_strain"><a class="viewcode-back" href="../../crack.html#quippy.crack.ConstantStrainRate.apply_strain">[docs]</a>    <span class="k">def</span> <span class="nf">apply_strain</span><span class="p">(</span><span class="bp">self</span><span class="p">,</span> <span class="n">atoms</span><span class="p">,</span> <span class="n">rigid_constraints</span><span class="o">=</span><span class="kc">False</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Applies a constant strain to the system.</span>

<span class="sd">        Parameters</span>
<span class="sd">        ----------</span>
<span class="sd">        atoms : ASE.atoms or quippy.Atoms</span>
<span class="sd">            Atomic configuration.</span>
<span class="sd">        rigid_constraints : boolean</span>
<span class="sd">            Apply (or not apply) strain to every atom.</span>
<span class="sd">            i.e. allow constrainted atoms to move during strain application</span>
<span class="sd">        &quot;&quot;&quot;</span>

        <span class="k">if</span> <span class="n">rigid_constraints</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span>
            <span class="n">initial_constraints</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">constraints</span>
            <span class="n">atoms</span><span class="o">.</span><span class="n">constraints</span> <span class="o">=</span> <span class="kc">None</span>

        <span class="n">newpos</span> <span class="o">=</span> <span class="n">atoms</span><span class="o">.</span><span class="n">get_positions</span><span class="p">()</span>
        <span class="bp">self</span><span class="o">.</span><span class="n">adjust_positions</span><span class="p">(</span><span class="n">newpos</span><span class="p">)</span>
        <span class="n">atoms</span><span class="o">.</span><span class="n">set_positions</span><span class="p">(</span><span class="n">newpos</span><span class="p">)</span>

        <span class="k">if</span> <span class="n">rigid_constraints</span> <span class="o">==</span> <span class="kc">False</span><span class="p">:</span>
            <span class="n">atoms</span><span class="o">.</span><span class="n">constraints</span> <span class="o">=</span> <span class="n">initial_constraints</span></div></div>



</pre></div>

           </div>
           
          </div>
          <footer>
  

  <hr/>

  <div role="contentinfo">
    <p>
        &copy; Copyright 2008-2016, James Kermode

    </p>
  </div>
  Built with <a href="http://sphinx-doc.org/">Sphinx</a> using a <a href="https://github.com/rtfd/sphinx_rtd_theme">theme</a> provided by <a href="https://readthedocs.org">Read the Docs</a>. 

</footer>

        </div>
      </div>

    </section>

  </div>
  


  <script type="text/javascript">
      jQuery(function () {
          SphinxRtdTheme.Navigation.enable(true);
      });
  </script>

  
  
    
   

</body>
</html>